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Abstract 


In  this  study,  the  subsurface  storage  and  transport  of  a  Dense  Non- Aqueous  Phase 
Liquid  (DNAPL),  trichloroethylene  (TCE),  was  evaluated  using  a  numerical  model. 
DNAPLs  are  organic  liquids  comprised  of  slightly  water-soluble  chemicals  or  chemical 
mixtures  that  have  a  density  greater  than  water.  Many  DNAPLs,  such  as  TCE,  are  used 
as  solvents  by  the  DoD  and  industry.  The  improper  disposal  and  handling  of  these 
chemicals  has  led  to  long  term  contamination  of  groundwater.  In  the  subsurface, 
DNAPLs  may  pool  atop  low  permeability  layers,  and  even  with  the  removal  or 
destruction  of  most  DNAPL  mass,  small  amounts  of  remaining  DNAPL  which  have  been 
transported  into  the  low  permeability  layer  can  dissolve  into  flowing  groundwater  and 
continue  as  a  contamination  source  for  decades.  Recently  developed  models  assume  that 
transport  in  the  low  permeability  zones  is  strictly  diffusive;  however  field  observations 
suggest  that  more  mass  is  stored  in  the  low  permeability  zones  than  can  be  explained  by 
diffusion  alone.  This  mass  may  be  in  the  form  of  separate  phase  DNAPL  or  dissolved 
phase  chlorinated  aliphatic  hydrocarbon  (CAH).  One  explanation  for  these  field 
observations  is  that  there  is  enhanced  transport  of  dissolved  CAHs  and/or  DNAPL  into 
the  low  permeability  layers  due  to  cracking.  Cracks  may  allow  for  advective-dispersive 
flow  of  water  contaminated  with  dissolved  CAHs  into  the  layer  as  well  as  possible 
movement  of  pure  phase  DNAPL  into  the  layer.  In  this  study,  a  numerical  flow  and 
transport  model  is  employed  using  a  dual  domain  construct  (high  and  low  permeability 
layers)  to  investigate  the  impact  of  cracking  on  DNAPL  and  CAH  movement.  Using 
literature  values,  crack  geometry  and  spacing  were  varied  to  model  and  compare  three 
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scenarios:  (1)  CAH  diffusion  into  an  uncracked  low  permeability  clay  layer;  (2)  CAH 
advection-dispersion  into  cracks,  and  (3)  separate  phase  DNAPL  movement  into  the 
cracks.  For  each  scenario,  model  simulations  are  used  to  show  the  evolution  and 
persistence  of  groundwater  contamination  down  gradient  of  the  DNAPL  source  caused 
by  back  diffusion  of  the  contaminant  out  of  the  low  permeability  layer  into  flowing 
groundwater.  This  study  found  cracking  will  cause  an  increase  in  transport  and  storage 
of  TCE  in  low  permeability  layers,  resulting  in  down  gradient  concentrations  above 
levels  of  concern  for  decades.  Further,  DNAPL  phase  TCE  within  cracks  can 
significantly  contribute  to  down  gradient  concentrations;  however,  the  extent  of  this 
contribution  is  very  dependent  upon  the  rate  of  DNAPL  dissolution.  Given  these 
findings,  remediation  goals  may  be  difficult  to  meet  if  source  remediation  strategies  are 
used  which  do  not  account  for  the  effect  of  cracking  upon  contaminant  transport  and 
storage  in  low  permeability  layers. 
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MODELING  THE  IMPACT  OF  CRACKING  IN  LOW  PERMEABILITY 

LAYERS  IN  A  GROUNDWATER  CONTAMINATION  SOURCE  ZONE  ON 
DISSOLVED  CONTAMINANT  FATE  AND  TRANSPORT 

1.  Introduction 

1.1  Background 

Widespread  use  of  chlorinated  solvents  such  as  tetrachloroethene  (PCE)  and 
trichloroethene  (TCE)  in  industrial  operations  over  the  last  century  has  resulted  in 
extensive  groundwater  contamination.  Poor  handling  and  disposal  of  these  chlorinated 
solvents  has  led  to  a  multitude  of  contaminated  sites  throughout  both  the  DoD  and 
industry.  The  EPA  estimates  that  over  60%  of  Superfund  sites  are  contaminated  with 
chlorinated  solvents.  If  these  sites  are  in  close  proximity  to  a  water  supply,  those  who 
ingest  the  water  are  at  an  increased  risk  for  developing  liver  problems  and  cancer  (EPA, 
2011).  The  solubility  of  many  chlorinated  solvents  may  be  as  high  as  several  g/L,  which 
exceeds  the  drinking  water  standard  by  a  factor  of  106  (EPA,  2011).  The  Environmental 
Protection  Agency  has  established  a  maximum  contaminant  level  (MCL)  for  TCE  in 
drinking  water,  the  contaminant  that  is  the  focus  of  this  study,  of  5  parts  per  billion,  or 
5pg/L. 

When  chlorinated  solvents  are  spilled  or  leaked  onto  the  ground,  they  move 
downward  as  a  dense  separate  phase  immiscible  liquid  or  “DNAPL”  (dense  nonaqueous 
phase  liquid).  Since  DNAPLs  have  a  specific  gravity  greater  than  water,  when  they  reach 
the  water  table,  they  continue  their  downward  migration  until  they  encounter  low 
permeability  layers.  The  DNAPL  will  form  pools  atop  these  layers.  These  pools  serve  as 
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a  persistent  contaminant  source  as  the  chlorinated  compound  dissolves  into  the  flowing 
groundwater.  In  this  study,  we  will  refer  to  dissolved  phase  chlorinated  compound  as  a 
“CAH”  (chlorinated  aliphatic  hydrocarbon).  Also,  as  the  DNAPL  migrates  downward, 
small  amounts  of  residual  DNAPL  are  left  behind  in  the  pore  spaces  between  the  aquifer 
solids.  This  residual  DNAPL  serves  as  another  source  of  persistent  contamination  as  it 
dissolves  into  groundwater.  The  flowing  groundwater  which  transports  the  CAH  will 
form  a  plume.  Plumes  vary  in  length  depending  upon  the  total  mass  of  the  contaminant, 
contaminant  properties  and  aquifer  conditions.  Water  will  flow  quickly  through  the  high 
permeability  layers  and  slowly  through  the  low  permeability  layers.  Figure  1  is  a 
conceptual  model  showing  DNAPL  distribution  in  the  subsurface,  as  well  as  the  plume 
that  forms  as  groundwater  flows  past  the  DNAPL  pools  and  residual. 


dissolved  and  sorbed 
contaminant 

bedrock 


by  Jeff  Heiderscheidt 


Water  Table 


residual  DNAPL 


plume 


Confining  Layer 

Low  Permeability  Layer 


Ground 
Water  Flow 
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Figure  1-A:  DNAPL  Distribution  in  the  Subsurface  (after  Heiderschiedt,  2010) 
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The  DNAPL  will  tend  to  migrate  down  through  the  high  permeability  layers,  pool 
atop  the  low  permeability  layers,  and  slowly  enter  the  low  permeability  layers  either  as  a 
CAH  or  DNAPL.  During  source  zone  remediation,  various  technologies  may  be  used  to 
remove  or  destroy  the  DNAPL  pools  and  residual.  However,  even  if  all  the  separate 
phase  residual  and  pooled  DNAPL  is  removed  from  the  high  permeability  zones,  if 
enough  contaminant  is  stored  within  the  low  permeability  layer  it  can  continue  to  act  as  a 
long-term  contaminant  source.  The  contaminant  can  be  stored  in  the  low  permeability 
layer  in  the  dissolved  phase  or  the  DNAPL  phase.  These  long-term  sources  can  continue 
to  contaminate  drinking  water  as  well  as  extend  the  cost  and  timeline  to  achieve 
remediation. 

As  described  above,  DNAPL  and/or  dissolved  chlorinated  compound  stored  in  low 
permeability  lenses  and  layers  in  the  subsurface  can  create  a  persistent  source  of 
contamination.  The  cleanup  of  low  permeability  lenses  is  very  difficult  and  often  long¬ 
term  contamination  continues  to  exist  at  sites  after  the  remediation  is  considered 
complete  (AFCEE,  2007).  The  movement  into  these  low  permeability  layers  is  typically 
modeled  as  Fickian  diffusion,  with  the  diffusion  coefficient  modified  to  account  for  the 
tortuosity  of  the  low  permeability  material,  as  well  as  for  retardation  due  to  contaminant 
sorption  to  the  solids  making  up  the  layer  (Parker  et  al.,  2008).  Recently,  however,  the 
applicability  of  this  Fickian  model  has  been  questioned,  and  the  potential  for  enhanced 
transport  in  these  low  permeability  layers  is  being  studied  (Miniter,  2011). 

One  hypothesis  for  the  enhanced  transport  is  that  the  low  permeability  materials 
contain  cracks.  Cracks  naturally  occur  in  clay  layers  and  can  be  caused  by  releases  of 
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pressure  due  to  erosion,  excavation,  or  changes  in  water  table  levels  (McKay  et  al., 

1993).  Cracking  is  dependent  upon  site  formation,  lithology,  and  composition.  Another 
hypothesis  is  that  the  interaction  of  the  DNAPL  mixture  with  low  permeability  lenses  in 
the  contamination  source  areas  can  result  in  an  alteration  of  the  properties  and  physical 
structure  of  the  low  permeability  lenses  (Demond,  2010).  This  alteration  may  also  lead  to 
cracking.  Thus,  cracking  of  low  permeability  materials,  whether  due  to  natural  processes 
or  the  interaction  of  a  DNAPL  mixture  with  the  low  permeability  lenses,  may  result  in 
enhanced  transport  of  contaminant  into  (and  out  of)  the  lenses.  This  enhanced  transport 
may  result  in:  (1)  advective  transport  of  dissolved  solvent,  (2)  DNAPL  entry  into  the 
cracks,  and/or  (3)  enhanced  diffusion  of  dissolved  solvent  into  the  cracks,  as  the  cracks 
have  lower  tortuosity  than  the  surrounding  matrix. 

It  has  been  reported  when  a  NAPL  is  in  contact  with  low  permeability  clay  layers  that 
the  hydraulic  conductivity  of  the  layers  can  increase  by  one  to  five  orders  of  magnitude. 
This  increase  in  hydraulic  conductivity  has  been  ascribed  to  interlayer  compression 
(Brown  and  Thomas,  1987).  The  shrinking  of  the  clay  layers  leads  to  the  formation  of 
cracks  and  micro  fractures  and  a  concomitant  increase  in  hydraulic  conductivity.  If  either 
dissolved  or  pure  phase  DNAPL  enters  these  cracks,  diffusion  into  the  low  permeability 
matrix  will  greatly  increase  due  to  a  larger  contact  area.  Both  naturally  occurring  cracks 
and  cracks  that  are  the  result  of  DNAPL  interaction  can  be  classified  by  aperture  size, 
depth,  surface  geometry,  surface  markings,  fabric  classification,  and  spacing  (Denness 
and  Fookes,  1969).  This  work  will  only  include  the  effects  of  crack  aperture,  depth,  and 
spacing. 
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In  order  to  quantify  the  impact  of  cracking  on  contaminant  transport,  a  dual-domain 
model  was  developed.  The  model  represented  cracks  as  a  so-called  “mobile  domain”, 
with  transport  of  dissolved  DNAPL  controlled  by  advection,  dispersion,  and  diffusion. 
The  clay  matrix  was  represented  as  an  immobile  domain,  where  diffusion  and 
equilibrium  sorption  controlled  transport  (Miniter,  2011).  A  first-order  rate  constant 
described  dissolved  DNAPL  transport  between  the  two  domains.  A  model  scenario  was 
constructed  where  a  pool  of  DNAPL  sat  within  a  high  permeability  sand  layer  atop 
cracked  clay  for  a  period  of  time.  The  model  was  used  to  simulate  concentrations  as  a 
function  of  time  down  gradient  of  the  DNAPL  source.  It  was  shown  that  the  existence  of 
cracks  in  the  clay  led  to  increased  concentrations  of  dissolved  DNAPL  downgradient, 
well  after  the  source  had  been  removed  (Miniter,  2011). 

1.2  Research  Objective 

The  primary  objective  of  the  research  is  to  model  the  impact  of  cracks,  either 
naturally  occurring  or  due  to  the  interaction  between  DNAPLs  and  low  permeability 
lenses  in  contamination  source  zones,  focusing  on  the  storage  and  transport  of  chlorinated 
solvents  within  these  lenses,  and  the  subsequent  impact  on  downgradient  dissolved 
contaminant  concentrations.  A  model  that  simulates  enhanced  diffusion  into  low 
permeability  lenses  that  was  previously  developed  by  Miniter  (2011)  will  be  further 
developed  to  model:  (1)  diffusion  only  into  the  cracks  and  surrounding  matrix,  and  (2) 
separate  phase  DNAPL  transport  into  the  cracks.  Results  from  these  simulations  will  be 
compared  with  earlier  conceptualizations  that  assume:  (1)  diffusion  into  uncracked  clay 
(AFCEE,  2007),  and  (2)  advective/dispersive  transport  in  cracks  and  diffusion  in  clay 
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(Miniter,  2011).  Properties  of  the  cracks  and  matrix  will  be  found  in  the  literature  and 
incorporated  in  the  model.  The  model  will  be  applied  to  simulate  changes  in  dissolved 
plume  behavior  resulting  from  cracks,  that  may  either  be  naturally  occurring  or  DNAPL 
induced,  in  low  permeability  clay. 

1.3  Research  Questions 

1 .  What  are  the  typical  characteristics  of  existing  cracks  in  low  permeability 
layers? 

2.  What  is  the  effect  of  cracking  on  the  transport  of  contaminants  into  and  out  of 
low  permeability  layers? 

3.  What  mathematical  models  can  be  used  to  simulate  transport  (e.g.  advection, 
dispersion,  diffusion,  and  DNAPL  transport)  into  and  out  of  cracks  in  low 
permeability  layers? 

4.  Compared  to  an  uncracked  source  zone,  how  is  the  flow  different  at  a  cracked 
source  zone? 

5.  What  is  the  effect  of  enhanced  transport  into  low  permeability  layers  on 
dissolved  plume  longevity  and  evolution? 

1.4  Research  Methodology 

1 .  The  initial  phase  of  the  study  involved  a  literature  review  to  a)  determine  if 
significant  cracking  occurs  naturally  in  low  permeability  layers,  b)  obtain 
parameters  to  characterize  these  naturally  occurring  cracks,  and  c)  determine  the 
appropriate  model  to  integrate  these  parameters  into  flow  and  transport  equations. 
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2.  Expand  the  existing  Miniter  (201 1)  model  to  include  diffusion  and  pure  phase 
DNAPL  movement  into  cracks  as  mass  transport  processes. 

3.  Use  the  Miniter  (201 1)  and  AFCEE  (2007)  models,  as  well  as  the  extended 
models  developed  in  the  previous  step,  to  quantify  and  compare  the  effects  of 
cracking.  The  comparison  will  consider  the  following  transport  processes:  (1) 
diffusion  only  into  cracks,  (2)  advection-dispersion  into  cracks,  (3)  pure  phase 
DNAPL  movement  into  the  cracks,  and  (4)  diffusion  into  uncracked  clay. 

1.5  Scope  and  Limitations  of  Research 

The  modeling  done  in  this  research  examines  highly  idealized  systems  to  examine 
the  possible  impact  of  cracking  in  subsurface  systems.  This  research  should  not  be  used 
to  predict  specific  concentrations,  rather  it  provides  a  qualitative  understanding  of 
DNAPL  plume  behavior.  The  values  used  for  simulations  are  based  on  common  values 
and  trends  found  in  literature,  not  a  specific  site.  The  model  used  in  this  study  assumes 
(1)  no  degradation  or  sorption  of  the  contaminant,  (2)  the  subsurface  material  properties 
in  each  layer  are  homogeneous  with  respect  to  space  and  time,  (3)  steady  state  flow,  and 
(4)  that  cracks  can  be  effectively  simulated  with  slight  changes  to  properties  in  a 
homogeneous  medium.  These  assumptions  are  necessary  for  both  model  execution  and 
to  create  comparisons  between  different  scenarios.  The  breakthrough  curves  and  mass 
balance  analyses  presented  are  presumed  to  be  an  accurate  comparison  of  the  effects  on 
down  gradient  plume  concentration  under  different  cracking  scenarios. 
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1.6  Definitions 


Advection  -  Flow  as  the  result  of  an  externally  applied  pressure  difference  or  as  a  result 
of  gravity/density  changes 

Basal  Spacing  -  Spacing  between  adjacent  layers  of  a  crystalline  structure 

Chlorinated  Aliphatic  Hydrocarbon  (CAH)  -  The  term  used  for  the  dissolved  phase 
DNAPL 

Crack  -  An  opening  in  a  material  caused  by  an  applied  stress  that  allows  fluid  to  freely 
enter  the  material 

Dense  Non  Aqueous  Phase  Liquid  (DNAPL)  -  A  fluid  which  has  a  density  greater  than 
water  and  is  also  relatively  immiscible  in  water 

Diffusion  -  Movement  of  a  dissolved  solute  governed  by  a  concentration  gradient 
according  to  Fick’s  Law 

Dispersion  -  Spreading  of  mass  due  to  spatial  and  temporal  variations  of  velocity  in  a 
flow  field.  In  this  work  dispersion  is  modeled  as  a  Fickian  process  to  capture  the 
heterogeneity  of  actual  porous  media 

Dissolution  -  The  process  of  dissolving  (e.g.,  from  a  DNAPL  phase  to  a  dissolved  phase) 

Enhanced  Diffusion  -  The  increased  amount  of  diffusion  into  porous  medium  than  can  be 
predicted  by  models  conventional  diffusion  models  governed  by  Fick’s  Law 
Entry  Pressure  -  The  pressure  required  for  a  DNAPL  to  enter  a  crack 
Hydraulic  Conductivity  -  The  capacity  of  a  medium  to  transmit  water 
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Maximum  Contaminant  Level  (MCL)  -  Legal  limit  for  a  contaminant  in  drinking  water 
set  forth  by  the  Environmental  Protection  Agency 

Permeability  -  The  ease  with  which  fluid  can  move  through  porous  material 
Retardation  -  Process  by  which  the  velocity  of  the  contaminant  becomes  less  than  the 
velocity  of  the  water  due  to  sorption 

Sorption  -  The  binding  of  a  contaminant  to  a  porous  medium 

1.7  Definition  of  Units 

In  this  work,  when  defining  variables  in  equations  the  units  will  follow  the 
variable  in  brackets.  A  list  of  units  is  shown  below. 

F  -  Force 

L  -  Length 

M  -  Mass 

T  -  Time 
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2.  Literature  Review 


2.1  Overview 

Understanding  the  processes  which  govern  contaminant  transport  in  aquifers  is 
vitally  important  in  setting  and  achieving  remediation  goals.  The  goal  of  this  research  is 
to  achieve  a  deeper  understanding  of  why  concentrations  down  gradient  of  sources  may 
remain  in  excess  of  MCLs  for  decades.  Due  to  gravity,  a  DNAPL  will  move  downward 
through  the  saturated  zone  of  porous  media  until  the  DNAPL  encounters  low 
permeability  layers  where  it  spreads,  forming  pools.  The  pools  of  DNAPL  will  slowly 
dissolve  into  flowing  groundwater,  as  well  as  diffuse  into  the  low  permeability  matrix. 
These  sources  create  dilute  plumes  which  can  extend  for  miles.  The  diffusion  into  the 
low  permeability  matrix  also  means  contaminants  will  diffuse  back  into  the  flowing 
groundwater,  perpetuating  the  plume,  even  after  the  pool  has  been  removed  or  completely 
dissolved.  This  chapter  examines  mechanisms  by  which  the  contaminant  is  transported 
into  the  low  permeability  matrix.  The  contaminant  may  enter  the  matrix  through  three 
possible  routes:  (1)  diffusion  of  the  CAH  into  the  uncracked  matrix,  (2)  advection, 
dispersion,  and  diffusion  of  the  CAH  into  the  cracked  matrix,  and  finally,  (3)  advection 
of  the  DNAPL  into  cracks  combined  with  diffusion  of  the  CAH  into  the  surrounding 
matrix. 

2.2  Aquifer  Characteristics 

An  aquifer  is  commonly  accepted  to  be  a  very  heterogeneous  medium.  Aquifers 
can  be  made  out  of  cracked  rock,  cobbles,  gravel,  sand,  clay,  silt,  or  most  commonly  a 
combination  of  many  different  materials.  Water  will  flow  through  an  aquifer  based  on 
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pressure  gradients  and  hydraulic  conductivity.  The  water  flows  through  the  medium  in 
accordance  with  Darcy’s  Law,  moving  from  areas  of  high  hydraulic  head  to  areas  of  low 
head.  This  means  water  can  flow  horizontally,  transversely,  and  vertically.  This  thesis 
considers  flow  through  sand  and  clay  only.  Based  on  the  relative  hydraulic  conductivities 
of  sand  and  clay,  water  will  flow  quickly  through  sand  and  extremely  slowly  through 
clay.  Sand  is  considered  a  high  permeability  medium  while  clay  is  considered  a  very  low 
permeability  medium. 

2.2.1  Naturally  Occurring  Cracking 

Cracking  occurs  naturally  in  low  permeability  layers,  and  these  cracks  may  allow 
for  enhanced  flow  and  transport  of  contaminants.  The  hydraulic  conductivity  of  cracked 
clay  is  commonly  two  to  three  times  higher  than  uncracked  clay  (McKay  et  al.,  1993). 
Cracks  can  be  formed  at  extensive  depths  from  weathering  or  stress  relief,  explaining 
their  occurrence  in  most  glacial  till  (Mackay  et  al.,  2000). 

Important  crack  characteristics  include  aperture,  spacing,  and  depth.  Crack 
apertures  generally  cannot  be  measured  directly,  therefore  many  investigators 
approximate  crack  apertures  in  field  and  laboratory  studies  using  the  cubic  law.  The  cubic 
law  uses  hydraulic  data,  and  an  important  assumption  is  that  the  crack  walls  are  two 
smooth  parallel  plates,  to  estimate  crack  aperture  (Sims  et  al.,  1996).  The  cubic  law  is 
shown  below  in  Equation  2. 1 : 

Qf  =  W^-W^l  Equation  2.1 

Qf  [L3T_1]  is  the  flow  rate  through  the  crack,  p  [ML-3]  is  the  fluid  density,  g  is 
gravitational  acceleration  [LT~2],  p  [ML'T1]  is  the  dynamic  viscosity  of  the  fluid,  2b  [L] 
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is  the  aperture  width,  W  [L]  is  width  of  the  sample,  AH  [L]  is  the  head  drop  over  the 
length  of  the  crack,  and  AL  [L]  is  the  crack  length  (Sims  et  al.,  1996). 

McKay  et  al.  (1993)  conducted  an  extensive  field  scale  investigation  in  order  to 
reliably  estimate  the  magnitude  and  distribution  of  crack  characteristics  in  a  clay  deposit. 
Previously  calculated  aperture  estimates  were  based  on  an  average  of  hydraulic 
conductivity  measurements  for  a  small  number  of  measurements.  McKay  et  al.  (1993) 
selected  the  Laidlaw  site  in  Lambton  County,  Ontario  due  to  the  extensive  knowledge 
base  established  by  previous  studies.  The  cubic  law  was  used  to  determine  crack 
aperture.  Transport  of  aqueous  contaminants  at  the  Laidlaw  site  was  expected  to  be 
governed  by:  (1)  advection  through  the  cracks,  (2)  diffusion  into  the  matrix  pore  water, 
and  (3)  retardation  and  degradation  processes  (sorption,  precipitation,  biodegradation)  in 
the  cracks  and  surrounding  matrix  (McKay  et  al.,  1993).  All  of  these  processes  are 
highly  dependent  upon  cracks.  Interestingly,  90%  of  the  cracks  in  the  upper  3.5m  had 
apertures  less  than  21  pm,  although  apertures  did  range  from  1  to  43  pm  (McKay  et  al., 
1993).  The  aperture  results  from  the  Laidlaw  site  were  similar  to  other  sites,  suggesting 
these  values  are  typical  for  glacial  till. 

A  study  conducted  by  O’Hara  et  al.  (2000)  estimated  the  size  and  variability  of 
crack  apertures.  The  methods  used  were  (1)  conventional  hydraulic  tests,  (2)  immiscible- 
phase  fluid  entry,  (3)  and  channel  identification  using  diffusion  halos  along  cracks.  The 
laboratory  study  used  a  column  which  was  0.5  m  in  diameter  and  0.5  m  in  length.  The 
column  was  extracted  from  between  3.7  and  4.2  m  depth  in  a  surficial,  slightly 
weathered,  clay  deposit.  Flowing  water  and  DNAPL  phase  TCE  were  used  to  identify 
areas  of  channeled  flow  in  the  cracks.  Although  horizontal  cracks  did  exist,  below  2  m 
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nearly  all  cracks  were  found  to  be  vertical  (O'Hara  et  al.,  2000).  The  mean  hydraulic 
conductivity  in  the  cracked  clay  was  found  to  be  three  times  greater  than  in  the  uncracked 
clay,  and  the  average  aperture  range  was  found  to  be  between  8  and  1 1  pm  (O'Hara  et  al., 
2000).  The  study  found  that  a  low  average  hydraulic  conductivity  does  not  necessarily 
mean  low  contaminant  transport  in  cracked  clay. 

Sims  et  al.  (1996)  studied  samples  from  a  weathered  cracked  clay  till  deposit 
using  a  flexible  permeameter.  The  samples  used  were  from  the  Laidlaw  Environmental 
Services  hazardous  waste  disposal  site  located  near  Sarnia,  Ontario.  The  cracks  at  the 
site  are  generally  attributed  to  desiccation  (Sims  et  al.,  1996).  The  clay  samples  were 
gathered  between  4  and  5  meters  below  land  surface.  Sample  locations  were  identified 
by  the  crack  halos  shown  in  Figure  2. 1 . 


Figure  2-A:  Diffusion  halos  surrounding  cracks  (Sims  et  al.,  1996) 

The  crack  halos  appeared  stained  due  to  oxidation  of  matrix  material.  The  highly 
oxidized  halos  are  indicative  of  recently  flowing  groundwater  (Sims  et  al.,  1996).  Two 
colors  of  halos  were  observed  at  the  site,  black/brown  and  grey/green.  The  black/brown 
coloration  of  the  halos  was  determined  to  indicate  active  flowing  oxygenated  ground 
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water,  which  resulted  in  staining  due  to  manganese  and  iron  oxidation.  The  grey/green 
halos  indicated  a  lesser  amount  of  oxygenated  groundwater  flow.  It  was  determined  the 
grey/green  halos  were  associated  with  dead  end  cracks,  meaning  offshoots  from  larger 
cracks  that  would  not  allow  significant  amounts  of  groundwater  movement.  The 
grey/green  halos  and  associated  cracks  were  most  abundant  at  a  depth  of  4  meters  and 
were  selected  for  study  (Sims  et  al.,  1996). 

Samples  were  gathered  by  using  a  70mm  Shelby  tube  over  isolated  halos.  The 
tube  was  driven  40  cm  into  the  till.  The  sample  was  sheared  from  the  surrounding  clay  to 
minimize  causing  an  increase  in  cracks.  The  samples  were  stored  at  4°C  and  sealed  in 
beeswax  to  prevent  further  desiccation.  Flow  tests  were  conducted  and  crack  apertures 
were  back  calculated  using  Equation  2.1,  the  cubic  law  relationship.  The  crack  apertures 
were  found  to  range  between  0  and  5  pm.  McKay  et  al.  (1993)  found  crack  apertures  to 
range  between  1  and  43pm  at  the  same  site,  and  Sims  et  al.  (1996)  suggested  that  the 
discrepancies  between  their  lab  based  findings  and  those  obtained  in  the  field  were 
because  samples  gathered  for  lab  testing  were  not  representative  of  field  conditions  due 
to:  (1)  the  larger  apertures  observed  in  the  field  may  have  become  plugged  or  filled  when 
extracted  and  transported  to  the  lab,  (2)  the  small  cracks  used  in  the  lab  do  not  represent 
field  scale  processes  such  as  flow  channeling,  and  (3)  cracks  which  were  open  in  the  field 
may  have  closed  during  sampling  or  lab  preparation. 

A  summary  of  crack  apertures,  depths,  and  spacing  from  various  studies  is  shown 
below  in  Table  2. 1 .  The  depth  listed  is  the  deepest  measurement  point  and  does  not 
necessarily  indicate  the  depth  at  which  the  crack  ends. 
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Table  2-1:  Crack  characteristics  found  in  the  literature 


Reference 

Depth 

(m) 

Spacing 

(m) 

Aperture 

(pm) 

D'Astous  et 
al.  (1989) 

<4 

0.04-0.1 

26-32 

Day  (1977) 

<18 

0.05-0.15 

1-14 

Grisak 

(1980) 

<7 

0.04 

4 

Henry  et  al. 
(1986) 

<16 

0.4 

50 

Hinsby  et  al. 
(1996) 

2-2.5 

- 

1-120 

Keller  et  al. 
(1986) 

12-18 

<.15 

11 

McKay  et  al. 
(1993) (2) 

1. 7-3.2 

0.04-0.13 

9-43 

McKay  et  al. 
(1993) 

<5 

0.02-1.0 

<43 

Pankow  et  al. 
(1984, 1986) 

<4 

0.03 

150 

O’Hara  et  al. 
(2000) 

3. 7-4.2 

- 

5-17 

Rudolph  et 
al.  (1991) 

<20 

1.5 

30 

Sims  et  al. 
(1996) 

4-5 

- 

1-5 

Thompson 

(1990) 

40-50 

1.2-5 

140-210 

2.2.2  DNAPL  Caused  Cracking 

While  cracks  occur  naturally  in  nature,  pooled  DNAPL  also  can  alter  the  physical 
properties  of  clay.  Many  studies  have  examined  the  impact  of  organic  liquids  on  basal 
spacing.  Basal  spacing  is  the  spacing  between  adjacent  layers  of  a  crystalline  structure, 
in  general  as  basal  spacing  increases  hydraulic  conductivity  decreases.  Clay  minerals  are 
considered  crystalline  structures.  A  clayey  deposit  will  generally  contain  a  significant 
amount  of  clay  minerals,  for  example  the  clay  mineral  content  at  the  aquitard  in  Dover 
AFB,  DE  ranged  between  18  and  35%  (Ball  et  al.,  1997). 
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Many  studies  have  examined  the  impact  of  organic  solvents  on  basal  spacing. 

One  study  by  Brown  and  Thomas  (1987)  studied  the  mechanism  by  which  an  organic 
liquid  can  change  the  hydraulic  conductivity  of  a  clay  by  measuring  basal  spacing, 
electrophoretic  mobility,  zeta  potential,  flocculation,  and  volume  change.  The  clays  used 
were  illite,  smectite,  and  kaolinite.  The  organic  liquids  used  in  this  study  were  acetone 
and  ethanol.  Brown  and  Thomas  (1987)  found  that  the  hydraulic  conductivity  was 
significantly  greater  in  kaolinitic  mixtures  exposed  to  acetone.  Further,  solutions  that 
were  70%  or  more  acetone  caused  the  same  effects  as  pure  acetone.  The  illitic  mixtures 
exhibited  small  increases  in  hydraulic  conductivity.  Dilute  acetone  solutions  (2-5%) 
caused  significant  increases  in  basal  spacing,  which  should  lead  to  a  decreased  hydraulic 
conductivity.  Brown  and  Thomas  (1987)  also  found  that  all  clays  when  in  contact  with 
an  organic  solution  would  swell  to  varying  degrees.  Brown  and  Thomas  (1987) 
concluded  that  as  organic  liquids  displace  water  in  equilibrium  with  clay  soil,  the  soil  will 
shrink,  causing  cracks  to  form.  These  cracks  can  act  as  channels  for  liquids  to  flow.  This 
results  in  an  increase  in  hydraulic  conductivity. 

Ayral  et  al.  (2011)  found  evidence  that  a  spilled  DNAPL  can  cause  cracking  in  an 
otherwise  unfractured  aquitard.  Three  hypotheses  are  proposed  by  the  researchers  for  the 
formation  of  cracks,  (1)  organics  liquids  decrease  the  basal  spacing  compared  to  water, 
"desiccating"  the  clay  (2)  surfactants  in  the  waste  can  change  the  wettability  (water  wet  to 
organic  wet),  enhancing  the  transport  of  the  organic  liquid,  and  (3)  interparticle  changes 
(floculation)  (Personal  Communication,  2012).  Only  the  first  hypothesis  is  discussed  in 
this  literature  review.  Ayral  et  al.  (2011)  compared  the  effects  of  various  organics  on 
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basal  spacing  including  the  pure  organic  liquid  solvent  TCE  with  a  waste  sample  of  TCE. 
Waste  DNAPLs  are  typically  mixtures  and  will  include  other  compounds  such  as 
surfactants.  The  changes  in  basal  spacing  of  montmorillonite  in  contact  with  various 
organics  were  examined  in  this  study.  Tests  showed  the  basal  spacing  was  smaller  when 
saturated  with  pure  TCE  and  waste  TCE  in  comparison  to  when  saturated  with  water 
(Ayral  et  al.,  2011).  The  basal  spacing  of  montmorillonite  was  examined  when  saturated 
with  a  surfactant.  While  the  basal  spacing  increased  as  compared  to  water,  the  basal 
spacing  of  the  waste  TCE  remained  closer  to  the  basal  spacing  of  pure  TCE.  This 
suggests  the  spacing  is  dominated  by  the  solvent  matrix  as  opposed  to  the  presence  of 
surfactants.  Interestingly,  Aryal  et  al.  (2011)  found  that  trichlorinated  organic  solvents 
cause  a  smaller  increase  in  basal  spacing  than  monochlorinated  or  dichlorinated  organics. 

In  the  laboratory,  clays  were  left  in  contact  with  TCE  waste  for  an  extended 
period  of  time.  Figure  2.2  (a)  depicts  a  crack  forming  after  being  exposed  to  TCE  waste 
for  ten  days,  and  Figure  2.2  (b)  depicts  increased  cracking  after  being  exposed  to  TCE 
waste  for  fifty  days.  Figure  2.3  displays  the  vertical  growth  of  cracks  (Personal 
Communication,  2012). 
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Figure  2-B:  (a)  Crack  formation  at  10  days,  (b)  crack  formation  at  50  days 


.  J 

Figure  2-C:  Vertical  Crack  Formation 
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2.3  Diffusion  of  the  CAH 


2.3.1  Conceptual  Model 

Transport  of  CAHs  into  low  permeability  layers  is  assumed  to  be  governed  by 
simple  Fickian  diffusion,  meaning  CAH  molecules  will  move  preferentially  from  zones 
of  high  concentration  to  low  concentration  due  to  Brownian  motion.  Mass  transfer  into 
the  low  permeability  layers  will  occur  due  to  the  relatively  high  concentrations  of  CAH 
found  around  the  DNAPL  pool  and  residuals.  The  net  result  of  the  relative  immobility  of 
the  pore  water  in  the  low  permeability  matrix  is  long-term  contaminant  accumulation  and 
high  contaminant  mass  storage. 

Back  diffusion  has  been  hypothesized  to  cause  long-term  plume  persistence.  The 
term  back  diffusion  refers  to  diffusion  of  contaminants  out  of  low  permeability  layers 
into  adjacent  high  permeability  zones.  Due  to  lower  concentration  gradients  driving  back 
diffusion,  back  diffusion  is  thought  to  be  a  much  slower  process  than  the  initial  diffusion 
process  (Sale  et  al.,  2008).  Although  back  diffusion  is  slow,  it  can  contribute  enough 
mass  to  the  flowing  groundwater  to  propagate  plumes  in  excess  of  MCLs  for  decades. 

Parker  et  al.  (2008)  concluded  that  back  diffusion  from  one  or  a  small  number  of 
thin  clayey  layers  in  a  sand  aquifer  can  cause  down  gradient  concentrations  to  remain 
above  MCLs  for  many  years  after  source  containment  or  removal.  A  site  heavily 
contaminated  with  TCE  in  Florida  was  studied  after  a  site  remediation  failed  to  achieve 
results  predicted  by  calculations.  Parker  et  al.  (2008)  compared  three  hypotheses  for  the 
cause  of  plume  persistence  after  the  source  zone  had  been  removed.  These  three 
hypotheses  were:  (1)  incomplete  source  zone  removal,  (2)  DNAPL  occurrence  down 
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gradient,  and,  (3)  back  diffusion  from  one  or  more  thin  clay  layers.  In  their  study,  Parker 
et  al.  (2008)  eliminated  the  first  two  hypotheses,  leaving  back  diffusion  as  the  only 
plausible  hypothesis.  The  study  found  that  even  a  clay  layer  of  less  than  0.2  m  thickness 
can  cause  “plume  persistence  due  to  back  diffusion  for  several  years  or  even  decades 
after  the  flux  from  the  source  is  completely  isolated”  (Parker  et  al.,  2008).  The 
concentrations  observed  in  the  plumes  from  CAH  stored  in  these  thin  clay  layers  were 
found  to  exceed  the  published  MCLs. 


2.3.2  Analytical  Model 

To  investigate  the  movement  of  solutes  into  clay,  Johnson  et  al.  (1989)  examined 
cores  of  uncracked  clay  near  a  site  studied  by  Goodall  and  Quigley  (1977)  and  Crooks 
and  Quigley  (1984).  These  cores  were  exposed  to  contaminants  and  were  then  analyzed 
for  chlorides,  organics,  and  total  organic  carbon.  Concentration  profiles  for  each  core 
were  developed,  and  the  distribution  of  chloride  was  examined  to  analyze  impacts  of 
diffusion.  Chloride  concentrations  in  each  core  were  high  at  the  surface  where  the  clay 
was  in  contact  with  the  waste  and  decreased  as  clay  depth  increased.  These  profiles 
suggested  the  primary  mechanism  of  transport  was  diffusion  and  the  profiles  were 
modeled  using  Equation  2.2,  Fick’s  second  law, 

Yt  =  Deff  JY  Equation  2.2 

where  C  [ML'3]  is  concentration,  t  is  time,  Deff  [L2T'1]  is  the  effective  solute  diffusion 
coefficient,  and  z  is  the  vertical  direction.  Deff  is  given  by  Equation  2.3: 

Deff  =  D0x  Equation  2.3 
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The  free  solution  diffusion  coefficient  (Do)  cannot  be  accurately  used  in  diffusion 
equations  due  to  the  tortuous  path  particles  follow  in  porous  media.  The  tortuosity  (x) 
ranges  between  0  and  1  and  accounts  for  the  complex  and  indirect  flowpath  solutes  travel 
through  a  porous  medium.  While  the  chloride  reached  maximum  depths  of  83  cm, 
organic  solutes  were  not  detected  past  a  depth  of  15  cm.  Johnson  et  al.  (1989)  concluded 
the  differences  in  diffusion  between  the  chloride  and  organic  solutes  was  due  to  the 
sorption  of  the  organics  to  clay.  The  effect  of  sorption  is  related  to  both  soil  and 
chemical  properties.  Commonly  the  effects  of  sorption  to  the  immobile  clay  phase  can  be 
represented  by  a  retardation  factor  (Rimm)  [-]  given  by  Equation  2.4, 

Rimm  =  1  +  Equation  2.4 

0 

where  pb  [M’L"3]  is  the  clay  bulk  density,  0  [-]  is  the  porosity  of  clay,  and  IQ  [M_1L3]  is 
the  sorption  partition  coefficient.  Rimm  is  incorporated  into  Equation  2.4  to  give  Equation 
2.5  and  2.6. 


Deff  - 


D0t 


dC  _  D0t  d2C 


Equation  2.5 
Equation  2.6 


The  model  employed  by  Johnson  et  al.  (1989),  which  used  Equation  2.6,  produced  results 
consistent  with  observed  concentrations  in  uncracked  clay. 

Removing  the  remaining  DNAPL  from  the  source  at  a  contaminated  site  is 
commonly  the  first  step  in  site  remediation.  The  goal  of  source  cleanup  is  to  remove  the 
source  of  the  dissolved  contaminant  plume.  Sale  et  al.  (2008)  investigated  how  source 
removal  will  affect  down  gradient  plume  concentrations.  Laboratory  and  analytical 
modeling  experiments  were  conducted  with  a  source  that  had  a  DNAPL  accumulate 
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above  a  low  permeability  capillary  barrier  (Sale  et  al.,  2008).  A  laboratory  study  was 
conducted  using  two  common  DNAPLs:  PCE  and  TCE,  and  a  light  non-aqueous  phase 
liquid  (LNAPL),  methyl-tert-butyl-ether  (MTBE).  The  chemicals  were  introduced  into 
highly  idealized  laboratory  experimental  boxes  with  a  high  permeability  transmissive 
zone  above  a  low  permeability  clay  zone.  The  assumed  mechanism  of  transport  was 
diffusion  only.  The  sources  were  introduced  and  then  removed  after  25  days.  Mass 
transport  was  then  evaluated  for  another  58  days.  The  study  found  that  15-44%  of  the 
contaminant  mass  was  stored  in  the  low  permeability  clay  layer  (Sale  et  al.,  2008).  The 
analytical  model,  which  predicted  results  found  in  the  laboratory,  then  was  applied  to 
simulate  longer  periods  of  time  allowing  for  diffusion  out  of  the  low  permeability  clay 
layer.  In  the  analytical  modeling,  the  source  was  present  for  the  first  1000  days  and  then 
removed.  Down  gradient  concentrations  remained  high  for  over  20  years,  with  the 
concentration  in  the  media  lm  from  the  source  decreasing  much  more  rapidly  than  the 
concentration  100m  from  the  source  (Sale  et  al.,  2008).  This  suggests  the  plume  itself 
can  act  as  a  secondary  source  of  contamination,  causing  contamination  of  clay  down 
gradient  of  the  DNAPL  source  due  to  diffusion  from  the  dissolved  plume  into  the  clay. 

2.3.3  Numerical  Model 

Parker  et  al  (2004)  modeled  TCE  transport  in  an  unfractured  minimally  weathered 
silt  aquitard.  At  the  site,  TCE  DNAPL  accumulated  at  the  bottom  of  a  10m  thick  sand 
aquifer  atop  a  20m  thick  silt  aquitard.  An  understanding  of  TCE  transport  through  the 
aquitard  ws  needed  due  to  the  pumping  of  drinking  water  in  an  underlying  sand  aquifer. 

A  one  dimensional  numerical  diffusive  transport  model  was  applied  to  the  site  to 
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determine  the  maximum  depth  to  which  the  aqueous  TCE  fronts  reached.  The  hydraulic 
as  well  as  diffusive  properties  were  determined  using  a  one  year  in  situ  tracer  study.  The 
model  was  verified  examining  cores  at  the  site  that  had  been  exposed  for  35-45  years. 

The  cores  verified  the  uncracked,  unweathered  nature  of  the  aquitard  and  indicated  that 
transport  was  diffusion-dominated.  Two  scenarios  were  considered:  a  time  period  of 
1200  years  with  no  advection  through  the  aquitard  and  a  time  period  of  500  years  with  a 
strong  vertical  hydraulic  gradient  through  the  silt  aquitard.  The  model  predicted  the  TCE 
would  reach  the  underlying  aquifer,  however  the  concentrations  reaching  pumping  wells 
would  not  be  above  MCLs  (Parker  et  al.,  2004). 

Chapman  and  Parker  (2011)  investigated  the  ability  of  three  numerical  models 
(HydroGeoSphere,  FEFLOW,  and  MODFLOW/MT3DMS)  to  simulate  two  scenarios. 
The  two  scenarios  were  (1)  the  experimental  situation  from  Sale  et  al.  (2008)  discussed 
above  and  (2)  a  two  layer  system  with  an  aquifer  atop  an  aquitard  solved  with  the 
analytical  solution  from  the  same  study  (Chapman  and  Parker,  2011).  The  results  of  their 
study  indicated  numerical  models  can  capture  field  scale  diffusion  processes  given 
sufficient  site  data. 

2.4  Advection  and  Dispersion  of  the  CAH 
2.4.1  Conceptual  Model 

Studies  conducted  by  Goodall  and  Quigley  (1977)  and  Crooks  and  Quigley 
(1984)  attempted  to  fit  an  analytical  diffusion-based  model  to  data  gathered  from  cracked 
clay  underlying  a  municipal  landfill  in  Ontario  (Johnson  et  al.,  1989).  Both  studies 
resulted  in  unsatisfactory  fits  of  the  models  to  the  data.  The  poor  fits  were  attributed  to 
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the  dominant  advective  forces  through  the  cracks,  which  were  not  accounted  for  by  the 
models.  In  this  section,  models  are  presented  that  simulate  enhanced  transport  by 
accounting  for  flow  through  cracks  in  low  permeability  media.  Conceptually,  the 
transport  of  dissolved  contaminants  by  advection  through  cracks  results  in  contaminant 
being  distributed  and  stored  deeper  within  the  low  permeability  layer  than  if  diffusion 
was  the  only  relevant  transport  process  . 

2.4.2  Analytical  Model 

Rowe  and  Booker  (1990)  used  an  analytical  model  to  examine  transport  through  a 
cracked  clay  landfill  liner  from  the  Ontario  municipal  landfill  studied  by  Goodall  and 
Quigley  (1977)  and  Crooks  and  Quigley  (1984).  Previously  the  authors  had  determined 
in  reviewing  the  literature  that  one  of  the  major  assumptions  used  in  selecting  clay  liners 
to  contain  landfill  waste,  that  unweathered  till  was  uncracked,  was  invalid  (Rowe  and 
Booker,  1990).  Rowe  and  Booker  (1990)  found  evidence  that  unweathered  till  could  in 
fact  be  cracked  to  extensive  depths,  which  would  allow  for  contaminant  transport.  Rowe 
and  Booker  (1990)  developed  a  model  that  assumed  1-D  contaminant  transport  in  the 
cracks  and  2-D  diffusion  into  the  surrounding  porous  medium.  The  parameters 
considered  by  the  model  included  pool  height,  concentration,  crack  depth,  porosity,  crack 
spacing,  and  crack  aperture.  The  processes  modeled  included  diffusion,  retardation,  and 
hydrodynamic  dispersion.  The  crack  spacings  ranged  from  0.5  to  5  meters,  and  the 
apertures  ranged  from  4  to  9  pm.  The  results  of  the  analytical  model  focused  on 
quantifying  the  effect  of  crack  spacing  and  Rowe  and  Booker  (1990)  concluded  that  even 
if  the  bulk  hydraulic  conductivity  is  known,  crack  density  can  significantly  impact  arrival 
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time  to  an  underlying  aquifer.  This  result  suggests  that  bulk  hydraulic  conductivity  is  not 
a  good  measure  of  contaminant  transport.  Rowe  and  Booker  (1990)  also  found  that  the 
concentration  of  contaminant  reaching  the  aquifer  was  very  sensitive  to  porosity  and  the 
depth  of  the  cracks. 

Ciahn  and  Tyner  (2011)  developed  a  2-D  radial  analytical  solution  for  solute 
transport  within  a  macropore  matrix  system  assuming  transverse  diffusion  and  advection 
of  contaminants  in  high  permeability  matrices  to  be  the  primary  transport  mechanism-^ 
Their  research  focused  on  a  small  conceptual  model  of  field  scale  back  diffusion.  The 
goal  of  their  research  was  to  determine  analytical  solutions  for  three  boundary  conditions 
which  would  be  validated  based  on  previously  published  data.  These  were:  (1) 
instantaneous  release  of  solute  into  a  macropore,  (2)  a  constant  concentration  of  solute  at 
the  top  of  the  macropore,  and  (3)  a  pulse  release  of  solute  into  the  macropore  (Cihan  and 
Tyner,  2011).  The  solutions  generated  by  Cihan  and  Tyner  (201 1)  assume  solute 
transport  at  the  macropore  level  is  governed  exclusively  by  advection,  and  solute 
transport  within  the  matrix  is  governed  exclusively  by  radial  diffusion.  Cihan  and  Tyner 
(201 1)  began  with  Equations  2.7  and  2.8: 


Equation  2.7 


rm9mdr  r_r™ 


Equation  2.8 


and  through  application  of  equations  describing  boundary  conditions  (1),  (2),  and  (3) 
above,  produced  solutions.  The  analytical  solutions  were  then  compared  with 
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experimental  data.  The  solutions  predicted  the  experimental  data  well,  although  accuracy 
was  found  to  be  dependent  upon  column  length. 

2.4.3  Numerical  Model 

A  numerical  model  developed  by  Grisak  and  Pickens  (1980)  used  an  advective- 
dispersive  transport  model  to  simulate  dissolved  solute  transport  in  cracked  media. 

Grisak  and  Pickens  (1980)  analyzed  the  impact  of  different  parameters  on  concentration 
histories  and  concentration  profiles  of  solutes  in  a  crack.  Important  parameters  were 
aperture  size,  water  velocity,  porosity,  and  the  distribution  coefficient.  Their  conceptual 
model  included  mobile  water  in  high  permeability  cracks  and  immobile  water  in  the  low 
permeability  matrix  surrounding  the  crack.  The  model  included  three  major  principles, 
(1)  diffusion  of  the  solute  into  the  low  permeability  matrix  from  the  high  permeability 
crack  (2)  advection  and  dispersion  due  to  flow  in  the  crack  and  (3)  linear  equilibrium 
sorption  in  the  matrix. 

Models  developed  prior  to  the  Grisak  and  Pickens  (1980)  model  emphasized  the 
total  effect  of  cracks  on  the  effective  permeability,  ignoring  the  interaction  between  the 
high  and  low  permeability  matrixes.  The  dominant  transport  mechanism  in  the  low 
permeability  no  flow  matrix  is  diffusion  and  the  dominant  transport  mechanisms  in  the 
high  permeability  cracks  are  advection  and  dispersion.  The  crack  was  modeled  as  shown 
in  Figure  2.4. 
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Figure  2-D:  Conceptual  model  of  crack  used  in  simulations  (Grisak  and  Pickens,  1980) 

A  constant  source  Co  was  maintained  at  the  upper  boundary,  with  a  zero 
concentration  gradient  at  the  center  of  matrix  blocks.  Assuming  a  constant  concentration 
of  C=Co  at  the  upper  boundary  (x=0)  would  simulate  the  effect  of  having  pooled 
DNAPL  atop  a  low  permeability  lense.  The  hashed  area  shows  the  diffusion  of  the 
contaminant  into  the  low  permeability  matrix  surrounding  the  crack.  The  model  allows 
for  solute  to  back  diffuse  out  of  the  matrix  into  the  flowing  crack  depending  on  the 
concentration  gradient. 

The  effect  of  diffusion  within  the  matrix  was  quantified  by  plotting  breakthrough 
curves  of  the  relative  concentration  (C/Co)  at  x=-0.76  m  for  diffusion  coefficients  ranging 
between  0.0  cm2/s  to  10~6  cm2/s.  As  the  diffusion  coefficients  increased  the  effect  of 
matrix  diffusion  became  more  pronounced.  The  authors  theorized  that  if  the  upper 
boundary  condition  of  C=Co  was  replaced  with  a  condition  of  C=0  the  solute  would 
diffuse  out  of  the  matrix  into  the  crack.  To  quantify  the  impact  of  aperture  size  on  solute 
mass  transfer  a  constant  velocity  in  the  crack  was  assumed  for  a  variety  of  apertures 
(Grisak  and  Pickens,  1980).  Reducing  the  aperture  size  reduces  the  quantity  of  solute 
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transported  in  the  crack  and  increases  the  relative  amount  of  solute  that  enters  the  matrix, 
since  for  a  given  diffusion  coefficient,  the  mass  flux  of  solute  into  the  matrix  is  controlled 
by  the  concentration  gradient  only.  When  the  aperture  size  is  decreased  the  fraction  of 
solute  transported  within  the  crack  is  decreased  while  the  fraction  of  solute  diffused  into 
the  matrix  is  increased  (Grisak  and  Pickens,  1980).  Increasing  the  velocity  in  the  crack 
increases  the  solute  flux.  Higher  velocities  produced  earlier  breakthrough;  however 
matrix  diffusion  became  a  significant  factor  after  a  short  time  in  all  scenarios  (Grisak  and 
Pickens,  1980).  Varying  the  matrix  porosity  produced  the  intuitive  result  that  larger 
matrix  porosities  resulted  in  greater  solute  transport  into  the  matrix  (Grisak  and  Pickens, 
1980).  The  sorption  coefficient,  IQ,  which  quantifies  the  relationship  between  dissolved 
and  sorbed  contaminant  concentrations  was  varied  in  the  matrix  only.  The  larger  the  IQ 
the  greater  the  solute  flux  into  the  matrix  (Grisak  and  Pickens,  1980). 

2.5  Advection  of  the  DNAPL  into  Cracks 

As  will  be  discussed  below,  evidence  from  both  the  laboratory  and  the  field,  as 
well  as  modeling  analyses,  suggest  that  DNAPL  sitting  atop  a  cracked  low  permeability 
layer  can  enter  the  cracks  when  the  entry  pressure  is  exceeded.  Diffusion  into  the  low 
permeability  matrix  will  then  occur,  as  the  DNAPL  dissolves  from  the  crack.  As  the 
DNAPL  mass  is  reduced  due  to  this  dissolution,  it  will  be  replenished  by  the  pool. 

2.5.1  Evidence  and  Conceptual  Model 

The  laboratory  experiment  conducted  by  O’Hara  et  al.  (2000)  which  examined 
fracture  flow  in  laboratory  columns  found  that  between  5  and  15%  of  cracks  contribute  to 


28 


DNAPL  flow,  and  all  other  cracks  could  contribute  to  CAH  advection  (O'Hara  et  al., 
2000). 

Hinsby  et  al.  (1996)  studied  a  clay  rich  till  deposit  located  near  Skaelskor, 
Denmark.  The  samples  were  taken  between  2  and  2.5  meters  below  land  surface.  The 
study  reports  crack  values  using  four  methods:  (1)  hydraulic  data,  (2)  solute  transport 
data,  (3)  colloid  transport  data,  and  (4)  measurements  of  nonwetting  fluid  entry  pressure 
for  a  DNAPL,  creosote.  Cracks  were  identified  by  grey/brown  halos  which  surrounded 
the  cracks,  and  a  sample  was  excavated  using  hand  tools  in  accordance  with  Jorgensen 
and  Spliid  (1992)  and  Jorgensen  and  Foged  (1994). 

In  method  (4),  creosote  was  added  to  the  soil  and  a  picture  of  a  sample  cross 
section  is  shown  below  in  Figure  2.5  with  a  cartoon  showing  the  magnification  of 
creosote  distribution  shown  in  Figure  2.6. 
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Figure  2-E:  Creosote  stains  in  crack  (Hinsby  et  al.  1996) 
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Figure  2-F:  Magnification  of  creosote  distribution  in  crack  (Hinsby  et  al.  1996) 


The  experiment  proved  that  under  pressure  highly  viscous  DNAPLs  will  enter  the  clay 
matrix.  Crack  apertures  were  calculated  for  method  (4)  using  Equation  2.13  which  is 
based  on  the  entry  pressure,  Pe,  required  for  creosote  to  enter  the  crack: 

2b  _  2o^os£>  Equation  2.13 

Where  a  [FL‘ 1  ]is  the  interfacial  tension  and  0  [°]  is  the  contact  angle.  Hinsby  et  al. 

(1996)  concluded  that  cracks  in  low  permeability  layers  may  permit  downward  migration 


31 


of  contaminants.  Crack  apertures  obtained  using  method  (1)  were  35  pm  and  56  pm  with 
crack  spacing’s  of  0.05  meters  and  0.20  meters,  respectively.  The  crack  aperture 
obtained  using  method  (2)  was  58  pm,  method  (3)  was  13  to  120  pm,  and  method  (4)  1  to 
94  pm. 

Hinsby  et  al  (1994)  also  ran  simulations  to  fit  chloride  tracer  data  from  the 
gathered  samples.  The  solute  model,  CRAFLUSH,  was  used  to  simulate  a  2-D  cross 
section  of  column.  CRAFLUSH  is  a  flow  and  transport  analytical  solution  which 
approximated  a  crack  as  two  parallel  plates  and  allows  for  diffusion  into  the  surrounding 
matrix.  The  hydraulic  conductivity  in  the  model  was  fixed  at  the  value  observed  in  the 
laboratory,  l.lxl O'6  m/sec.  Crack  spacing  was  varied  until  the  simulated  curves  fit  the 
observed  data.  Two  methods  were  used  for  curve  fitting.  The  first  used  the  cubic  law  to 
determine  crack  aperture,  fixing  crack  spacing  and  varying  vertical  hydraulic 
conductivity.  The  second  method  fixed  the  crack  spacing  and  varied  crack  aperture. 

Both  methods  yielded  approximately  the  same  results  for  crack  aperture,  58  pm.  The 
results  from  the  second  method  also  suggested  that  crack  spacing  had  little  impact  upon 
solute  transport.  In  this  experiment,  the  minimal  influence  of  crack  spacing  is  most  likely 
due  to  the  high  velocity  of  flow  in  the  cracks  and  the  short  length  of  the  column  used  in 
the  experiment.  In  decreased  flow  scenarios  the  simulations  were  more  sensitive  to  crack 
spacing  (Hinsby  et  al.,  1996). 

2.5.3  Numerical  Model 

A  model  of  contaminant  transport  in  cracked  aquifers  was  presented  by  both 
Mackay  and  Cherry  (1989),  and  Kueper  and  McWhorter  (1991).  These  models  suggested 
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DNAPL  will  enter  cracks  which  provide  the  least  resistance  to  flow  and  will  continue  to 
flow  into  the  crack  until  capillary  forces  impede  DNAPL  flow.  In  these  cracks  some  of 
the  DNAPL  will  form  pools  at  the  bottom  of  cracks,  leave  residual  along  the  flow  path, 
and  enter  the  surrounding  matrix  through  diffusion.  Both  studies  found  that  the  crack 
porosities,  void  space  of  the  crack  per  total  volume  of  cracked  media,  are  commonly 
several  orders  of  magnitude  less  than  most  granular  aquifers,  and  therefore  little  DNAPL 
can  be  stored  in  the  cracks.  However,  this  also  means  flow  in  cracks  can  allow  for 
widespread  migration  of  contaminants.  A  revised  model  was  proposed  in  Parker  et  al. 
(1994;  1997)  which  demonstrated  that  the  ability  of  the  matrix  to  store  the  dissolved 
DNAPL  exceeds  the  ability  of  the  crack  to  store  the  DNAPL.  Once  the  DNAPL  enters 
the  crack  it  is  surrounded  by  a  thin  layer  of  water,  resulting  in  a  large  DNAPL/water 
interfacial  area  relative  to  the  DNAPL  volume.  This  allows  for  a  large  quantity  of  the 
DNAPL  to  dissolve  even  though  many  DNAPLs  have  relatively  low  solubilities.  Once 
dissolved,  the  contaminant  will  move  based  on  concentration  gradients  potentially 
driving  the  contaminant  further  into  the  surrounding  matrix. 

The  model  proposed  by  Parker  et  al.  (1994)  for  immiscible  organic  liquids  in 
cracked  porous  media  included  the  effects  of  diffusion  on  the  persistence  of  organic 
liquids  in  the  cracks.  The  model  takes  into  account  the  very  high  DNAPL  surface  area  to 
volume  ratio,  which  allows  for  fast  diffusion  at  the  NAPL/water  interface.  The  surface 
area  of  NAPL  in  cracks  is  large  compared  to  the  surface  area  of  a  pool  atop  a  low 
permeability  layer.  The  layer  of  water  surrounding  the  DNAPL  is  assumed  to 
instantaneously  reach  solubility  and  a  chemical  concentration  gradient  will  be  established 
with  the  surrounding  matrix.  As  the  DNAPL  in  the  cracks  dissolves  and  diffuses  into  the 
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matrix,  it  is  replenished  by  the  DNAPL  pool  sitting  atop  the  low  permeability  layer. 
Eventually,  the  DNAPL  in  the  cracks  will  become  disconnected  ganglia  as  the  height  of 
the  DNAPL  pool  (and  therefore,  the  pressure)  decreases  and  the  effects  of  diffusion  are 
more  pronounced  (Parker  et  al.,  1994).  Mass  transfer  from  the  DNAPL  to  the  immobile 
pore  water  is  included  in  the  Parker  et  al.  (1994)  model.  The  model  accounts  for 
reduction  in  the  volume  of  immiscible  DNAPL,  thereby  reducing  the  maximum 
penetration  depth.  Linally  when  all  or  most  DNAPL  dissolves,  the  CAH  continues  to 
migrate  into  the  matrix  as  long  as  the  concentration  gradients  favor  this  migration.  When 
the  pool  is  no  longer  atop  the  low  permeability  layer  and  clean  water  is  passing  over  the 
cracks  contaminants  will  be  removed  from  the  matrix  by  back  diffusion.  CAH  removal 
will  be  controlled  by  diffusion  and  desorption  from  the  matrix  solids  and  pore  water.  The 
time  required  for  the  above  processes  is  dependent  upon  the  properties  of  the  DNAPL, 
geologic  material,  and  quantity  of  DNAPL  (Parker  et  al.,  1994). 

Parker  et  al.  (1994)  also  studied  the  impact  of  crack  spacing  on  contaminant  flux. 
Previous  models  had  assumed  a  single  crack  in  an  infinite  porous  medium  which  allows 
large  concentration  gradients  to  form.  In  reality  concentration  gradients  will  form  around 
all  cracks,  and  when  the  concentration  gradients  compete,  the  diffusion  and  mass  flux 
from  crack  surfaces  will  decrease.  The  time  when  diffusion  profiles  will  meet  between 
cracks  after  a  release  of  DNAPL  is  dependent  upon  the  distance  between  cracks,  DNAPL 
properties,  and  matrix  properties.  Medium  to  narrow  aperture  size,  high  aqueous 
solubility,  large  porosity,  and  high  sorption  capacity  of  the  matrix  were  shown  to  enhance 
the  rate  of  the  NAPL  dissolution  (Parker  et  al.,  1994). 
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Kueper  and  McWhorter  (1991)  examined  conditions  for  DNAPL  transport  into 
cracks.  In  their  equations  it  is  assumed  the  matrix  is  initially  fully  saturated  with  water, 
and  that  water  is  the  wetting  fluid  and  the  DNAPL  is  the  nonwetting  fluid. 

Capillary  pressure  (Pc)  is  defined  as  the  difference  in  pressure  between  that  in  the 
nonwetting  fluid  (Pnw)  and  that  in  the  wetting  fluid  (Pw)  as  shown  in  Equation  2.9. 

Pc  =  Pnw  +  Pw  Equation  2.9 

As  previously  stated,  for  the  DNAPL  to  enter  the  clay  matrix  the  capillary  pressure  at  the 
top  of  the  crack  must  exceed  the  entry  pressure  (Pe)  of  the  DNAPL.  Two  equations  will 
predict  Pe,  Equation  2.10  is  for  irregular  crack  patterns  approximated  by  two  parallel 
plates,  and  Equation  2.11  is  for  circular  cracks. 

Pe  =  parallell  plates  Equation  2. 10 

Pe  =  circular  Equation  2.11 

a  [FD1]  is  the  interfacial  tension  between  the  DNAPL  and  the  water,  0  [-]  is  the  contact 
angle  measured  through  the  wetting  phase,  and  2b  [L]  is  the  crack  aperture.  From  the 
equations  it  is  noted  that  DNAPL  entry  pressure  is  inversely  proportional  to  crack 
aperture;  that  is,  larger  crack  apertures  require  lower  Pe  for  DNAPL  entry  into  the  crack. 
In  reality  the  entry  pressure  will  lie  in  between  the  entry  pressure  predicted  by  Equations 
2.10  and  2.1 1.  Figure  2.7  shows  two  possible  idealizations  of  a  rough  walled  crack. 
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Figure  2-G:  Rough  crack  and  two  idealizations  (Kueper  and  McWhorter,  1991) 

Cracks  in  nature  are  generally  elongated  and  irregular,  therefore  Equation  2.9,  which 
results  in  a  lower  entry  pressure,  will  be  used  for  calculations.  Using  a  lower  entry 
pressure  is  conservative.  The  DNAPL  pool  height  (HD)  required  to  create  a  given  entry 
pressure  is  given  by  equation  2.12: 

HD  =  —  Equation  2.12 

u  Apg2b  H 

'J 

where  Ap  [ML-  ]  is  the  density  difference  between  water  and  the  DNAPL.  Equation  2.12 
displays  an  inverse  relationship  between  aperture  width  and  the  pool  height  required  for 
DNAPL  entry.  The  Kueper  and  McWhorter  (1991)  model  also  indicated  that  DNAPL 
transport  through  a  cracked  aquitard  increased  with  downward  water  gradients  and 
decreased  with  upward  water  gradients  across  the  aquitard.  (Kueper  and  McWhorter, 
1991). 
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Reynolds  and  Kueper  (2002)  examined  the  effects  of  crack  aperture,  matrix 
porosity,  and  matrix  organic  carbon  on  the  migration  of  five  DNAPLs  using  numerical 
simulations.  The  study  found  that  aperture  is  the  most  important  factor  impacting 
migration  of  DNAPLs  in  cracks.  Particularly  in  large  cracks  matrix  diffusion  does  not 
retard  the  rate  of  DNAPL  migration.  By  increasing  the  crack  aperture  from  15  pm  to  50 
pm  there  was  a  20  fold  increase  in  the  rate  of  DNAPL  migration  downward  (Reynolds 
and  Kueper,  2002). 

Murphy  and  Thomson  (1993)  developed  a  dynamic  two-dimensional  two-phase 
flow  model  for  a  single  aperture.  The  system  was  a  finite  volume  implementation  of  the 
cubic  rule,  Equation  2.1,  and  assumes  incompressible  flow  between  parallel  plates. 
Esposito  and  Thomson  (1999)  extended  the  two  phase  crack  flow  model  developed  by 
Murphy  and  Thomson  (1993).  Their  numerical  model  includes  transient  two-phase  flow, 
non-equilibrium  dissolution,  advective-dispersive  transport  in  the  crack,  and  three- 
dimensional  matrix  diffusion  (Esposito  and  Thomson,  1999).  The  authors  were 
investigating  (1)  the  role  dissolution  and  diffusion  have  on  DNAPL  disappearance  and 
(2)  how  changing  flushing  rates  affect  DNAPL  mass  removal.  The  model  approximated 
two  phase  flow  in  a  single  crack  as  parallel  plate  flow.  Apertures  are  generally  rough  and 
can  vary  in  width  considerably,  therefore  Esposito  and  Thomson  (1999)  created  log 
normal  aperture  distributions  using  an  algorithm  developed  by  Robin  (1991).  This  allows 
the  domain  for  the  model  to  be  a  crack  network  with  varying  apertures  incorporating  the 
narrowing  of  cracks  with  depth.  Figure  2.8  is  an  example  of  an  Esposito  model  domain 
(2D  cross  section). 
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Figure  2-H:  Variable  aperture  model  domain  (Esposito  and  Thomson,  1999) 


Figure  2.8  has  an  average  aperture  value  of  210  pm  and  was  used  to  investigate  the  role 
of  dissolution  and  diffusion  on  DNAPL  disappearance.  As  seen  above  the  domain  is 
relatively  small  and  therefore  it  was  necessary  to  assume  that  the  size  of  the  control 
volume  (shown  above)  sufficiently  accounted  for  aperture  variations  (Esposito  and 
Thomson,  1999).  Also,  the  model  does  not  account  for  large  scale  heterogeneities  in  the 
media. 

The  authors  used  three  materials  for  their  simulations,  Clay  1  which  had  a 
porosity  of  0.35  and  a  fraction  of  organic  carbon  (foc)  of  0.01,  Clay  2  which  had  a 
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porosity  of  0.55  and  a  foc  of  0.005,  and  shale/sandstone  which  had  a  porosity  of  0.10  and 
a  foc  of  0.002.  The  retardation  coefficient  of  Clayl  and  2  was  6.8  while  the 
Shale/Sandstone  was  7.0.  Clay  1  and  2  used  the  same  diffusion  coefficient  while  the 
Shale/Sandstone  used  a  much  lower  coefficient  (Esposito  and  Thomson,  1999). 

The  first  investigation  compared  the  time  it  took  for  the  DNAPL  TCE  to  travel 
0.13  and  0.43  meters  in  the  porous  media.  In  Clay  1,  it  took  64.8  years,  in  Clay  2,  22.0 
years,  and  in  the  Shale/Sandstone  it  took  486.1  years  (Esposito  and  Thomson,  1999). 
These  results  mean  a  high  porosity  and  low  foc  allow  for  faster  DNAPL  travel.  A  smaller 
porosity  and  diffusion  coefficient  in  the  Sandstone/Shale  impaired  downward  movement 
but  also  prevented  significant  mass  movement  into  the  surrounding  matrix.  These 
findings  were  consistent  with  the  findings  in  Parker  et  al.  (1994). 

The  second  investigation  analyzed  the  effectiveness  of  mass  removal  methods 
which  rely  on  hydraulic  gradients.  20g  of  the  DNAPL  was  allowed  to  enter  the  grid  and 
then  three  different  hydraulic  gradients  (0.4,  0.04,  and  0.004)  were  used  to  flush  the 
DNAPL.  The  model  assumed  mass  removal  began  immediately  following  the  DNAPL 
release,  which  is  very  optimistic  since  most  spills  are  not  detected  until  many  years  after 
the  spill  occurred.  The  time  required  to  remove  99%  of  the  mass  ranged  from  10  days  to 
several  hundred  years,  which  means  remediation  technologies  that  use  hydraulic 
gradients  (such  as  pump  and  treat)  may  be  severely  impaired  in  a  cracked  clay 
environment.  Further,  if  the  pool  of  DNAPL  remains  as  a  long  term  source  the  high 
concentration  gradients  surrounding  the  cracks  will  drive  more  contaminant  into  the 
surrounding  matrix  making  complete  mass  removal  even  more  difficult  (Esposito  and 
Thomson,  1999). 
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3.  Methodology 


3.1  Overview 

Cracking  in  clay  is  hypothesized  to  contribute  to  enhanced  diffusion  and  storage 
of  contaminants  in  low  permeability  layers.  Cracks  are  known  to  form  naturally  in  low 
permeability  layers  as  a  result  of  natural  loading  and  unloading  cycles  as  well  as 
dessicationdesiccation  (McKay  and  Fredericia,  1995)  (Esposito  and  Thomson,  1999). 
Recent  research  has  also  indicated  pooled  DNAPL  can  cause  cracks.  Pooled  DNAPL  is 
hypothesized  to  then  enter  the  cracks.  In  this  study,  a  numerical  grid  is  created  to 
simulate  a  high  conductivity  sand  layer  sitting  atop  a  low  conductivity  clay  layer.  The 
results  reported  in  Miniter  (201 1)  are  verified,  and  the  transport  of  TCE  is  modeled  in 
three  scenarios.  These  scenarios  are:  (1)  transport  into  uncracked  clay,  (2)  transport  of 
the  dissolved  TCE  into  cracked  clay,  and  (3)  transport  of  DNAPL  phase  TCE  into 
cracked  clay.  This  work  is  an  extension  of  the  work  done  by  Miniter  et  al.  (201 1)  and 
will  evaluate  the  impact  of  cracking  on  enhanced  diffusion  and  storage  into  low 
permeability  layers.  The  model  in  scenarios  (2)  and  (3)  assumes  the  existence  of  cracks 
in  the  low  permeability  layers,  but  does  not  differentiate  between  naturally  occurring  or 
DNAPL  caused  cracks. 


3.2  General  Description 

This  research  assumes  a  heterogeneous  system  composed  of  three  media:  sand, 
clay,  and  cracked  clay.  A  DNAPL  source  is  simulated  within  a  high  permeability  sand 
layer  atop  a  low  permeability  clay  layer.  The  contaminant  will  be  transported  into  the 
clay  through  various  transport  processes  including  diffusion,  advection  of  the  CAH,  and 
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advection  of  the  DNAPL.  In  uncracked  clay,  it  is  hypothesized  the  only  mechanism  of 
transport  of  the  DNAPL  into  the  clay  is  diffusion.  In  cracked  clay  it  is  hypothesized  the 
TCE  will  either  advect  as  a  pure  phase  DNAPL  or  as  dissolved  CAH  into  the  cracks.  In 
the  cracked  clay,  the  contaminant  will  then  diffuse  into  the  surrounding  clay  matrix  based 
on  concentration  gradients.  The  increase  in  transport  due  to  cracking  will  be  termed 
“enhanced  diffusion.”  The  enhanced  diffusion  into  the  clay  could  lead  to  a  significant 
increase  of  contaminant  storage  in  the  clay  layer.  Figure  3.1  is  a  conceptual  diagram 
depicting  the  source  zone. 


Dissolution  of  the  DNAPL  into  the  high  flow  sand  layer  will  remove  much  of  the 
DNAPL  mass  and  create  a  down  gradient  plume.  This  plume  will  result  in  concentration 
gradients  down  gradient  of  the  source  zone  which  will  permit  further  transport  by 
diffusion  into  the  low  permeability  clay  layer.  Once  the  source  is  removed,  it  is 
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hypothesized  the  CAH  will  diffuse  out  of  the  low  permeability  clay  layer  perpetuating  the 
plume  at  lower  concentrations.  Monitoring  wells  placed  in  the  high  permeability  layer 
will  show  the  concentration  history  of  the  dissolved  plume  down  gradient  of  the  source. 
Modeling  will  allow  us  to  quantify  the  mass  stored  in  the  low  permeability  layer,  as  well 
as  the  down  gradient  concentration  history.  Analysis  of  the  model  results  will  provide  an 
increased  understanding  of  the  relationship  between  cracking,  mass  storage,  and  plume 
behavior.  Figure  3.2  depicts  plume  formation  and  a  plume  observation  point. 


Figure  3-B:  Plume  Formation  and  Observation  Point  Conceptual  Diagram 
3.2.1  Assumptions 

Many  key  assumptions  are  required  to  model  contaminant  transport  in  an  aquifer 
system.  In  general,  the  subsurface  is  anisotropic  and  heterogeneous.  This  means  flow 
can  be  significantly  different  at  two  different  points  in  the  same  medium.  This  model 
assumes  the  same  hydraulic  and  chemical  properties  within  each  defined  layer.  The 
source  is  assumed  to  be  instantaneously  removed,  simulating  total  remediation  of 
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DNAPL  mass  in  the  high  permeability  zone.  Further,  when  modeling  the  cracks  it  is 
necessary  to  assume  uniform  aperture,  spacing,  and  depth.  When  modeling  DNAPL 
transport  into  the  cracks  it  was  assumed  the  DNAPL  has  fully  penetrated  the  crack  prior 
to  the  start  of  the  simulation.  Also,  to  ensure  the  DNAPL  transport  is  more  realistically 
modeled,  at  ten  years  once  the  source  is  removed  the  DNAPL  is  again  assumed  to  have 
fully  saturated  the  crack. 

In  this  thesis,  we  assume  cracked  clay  can  be  modeled  as  an  equivalent 
homogeneous  anisotropic  medium  with  increased  vertical  hydraulic  conductivity.  Figure 
3.3  below  provides  a  visual  representation  of  this  assumption. 


Figure  3-C:  Vertical  Flow  Approximation  for  Cracking 
3.3  Governing  Equations 

Many  of  the  equations  used  in  the  model  and  for  basic  calculations  are  the  same 
equations  used  by  Miniter  (2011).  For  a  further  explanation  of  equations  consult  the 
work  of  Minter  (2011). 
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3.3.1  Dissolved  Transport 

The  primary  equation  governing  contaminant  transport  in  porous  media  is  shown 
below  in  Equation  3.1: 

IrCE=^{Di,^y^ivlQ-pbdlt  Equal,  on3.1 

where  JTce  [ML'2T_1]  is  the  contaminant  flux,  x  [L]  is  the  distance  along  the  respective 

Cartesian  coordinate  axis,  C  [ML'  ]  is  the  concentration  of  dissolved  contaminant, ,  Djj 

2  1  1 

[L  T'  ]  is  the  hydrodynamic  dispersion  coefficient  tensor,  Vj  [LT‘  ]  is  the  linear  pore 

'J  'J 

water  velocity,  pb  [ML'  ]  is  the  bulk  density,  and  S  [ML'  ]  is  the  sorbed  contaminant. 
While  Equation  3.1  is  applied  throughout  the  entire  model  domain,  transport  in  sand, 
clay,  and  cracked  clay  will  vary  based  on  the  specific  medium’s  properties. 

The  hydrodynamic  dispersion  coefficient  combines  the  effects  of  mechanical 
dispersion  and  diffusion  into  one  term.  This  relationship  is  shown  in  Equation  3.2 
(Gupta,  2008): 

=  Deff  +  Dm  Equation  3.2 

Dm  is  the  mechanical  dispersion  coefficient  given  by  Equation  3.3: 

Dm  =  otjjVj  Equation  3.3 

Where  a,,  [L]  is  the  dispersivity  parameter  and  Vj  [LT1]  is  the  linear  pore  water  velocity. 
Deff  is  the  effective  molecular  diffusion  coefficient.  The  effective  diffusion  coefficient 
(Deff)  is  related  to  the  free-solution  diffusion  coefficient  (Do)  shown  below  in  Equation 
3.4. 


Deff  D0t 


Equation  3.4 
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where  x  is  a  tortuosity  factor  (0<x<14— that)  that  accounts  for  the  hindrance  to  diffusion 
through  porous  media  (Gupta,  2008). 

Diffusive  transport  in  cracks  and  into  the  surrounding  matrix  is  governed  by 
Fick’s  second  law  shown  below  in  Equation  3.5  (Johnson  et  al.,  1989). 

f  =  £>«//§  Equation  3.5 

Where  Deff  [L2T_1]  is  the  effective  diffusion  coefficient,  C  [ML'3]  is  the  aqueous 
concentration,  and  x;  [L]  is  the  respective  Cartesian  coordinate.  As  seen  in  Equation  3.5, 
the  transport  of  mass  into  the  clay  is  determined  by  the  concentration  gradients. 

TCE,  like  many  other  DNAPLs,  will  readily  adsorb  to  organic  material  in  porous 
media.  As  the  concentration  increases  in  the  flowing  groundwater,  the  mass  adsorbed  to 
organic  materials  will  increase  proportionally.  Sorption  may  be  modeled  as  a  kinetic 
process,  as  is  shown  below  in  Equation  3.6 

d  S 

pb  —  =  9a(C  —  kdS )  Equation  3.6 

Where  S  is  the  sorbed  mass  [M],  0  [-]  is  the  media  porosity,  kd  [MLwater'  ]  is  the  sorption 
coefficient,  and  a  [T1]  is  the  first  order  mass  transfer  coefficient  between  the  water  and 
solid  grains  of  the  aquifer.  In  this  study,  sorption  is  assumed  to  be  negligible  (a  =  0). 

The  DNAPL  within  the  cracks  will  slowly  diffuse  into  the  flowing  water  based  on 
concentration  gradients.  In  order  to  effectively  model  DNAPL  dissolution,  a  step 
function  must  be  used.  Equation  3.7  is  the  differential  equation  modeling  the  first  order 
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DNAPL  dissolution,  and  Equation  3.8  is  the  step  function  indicating  DNAPL  dissolution 
ends  when  the  DNAPL  is  completely  dissolved  and  no  longer  present. 


P^pi2£ffi  =  /?(C-f)  Equation  3.7 

F  =  f?  Equation  3.8 

Jnapl  — 

Where  Pnapl  [ML'3]  is  the  NAPL  density,  Snapl  [-]  is  the  NAPL  saturation,  Cs  [ML'3]  is 
the  DNAPL  solubility  in  water,  p  [T'1]  is  the  first  order  mass  transfer  coefficient  between 

'J 

the  DNAPL  and  the  flowing  water,  and  C  [ML'  ]  is  the  concentration  of  the  dissolved 
DNAPL.  P  is  determined  using  the  following  relationship  in  Equation  3.9  (Christ  et  al., 
2006): 


Pchrist 


Equation  3.9 


Where  ko’fT'1]  is  a  fitting  parameter,  M(t)  [M]  is  the  time  dependant  mass,  Mo  [M]  is  the 
mass  at  time  zero,  and  Pchrist  H  is  a  fitting  parameter.  For  this  work,  p  is  assumed  to  be 

constant,  and  its  value  determined  by  using  =  0.5.  Values  used  from  Christ  et  al. 

M0 

(2006)  for  ko’  and  Pchrist  are  8.2e-3  d'1  and  0.85,  respectively. 


3.3.2  Cracking 

The  number  of  cracks  in  the  clay  at  the  source  can  be  determined  using  Equation 

3.10: 

c  =  (^  +  1)(^  +  1)  Equation  3.10 
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where  c  [-]  is  the  number  of  cracks,  L  [L]  is  the  length  of  the  source  zone,  W  [L]  is  the 
width  of  the  source  zone,  and  2B  is  the  distance  between  cracks.  In  order  to  model 
advection  of  DNAPL  within  the  cracks  it  is  necessary  to  calculate  crack  volume.  The 
equation  for  a  single  crack  volume  is  shown  below  in  Equation  3.11: 

VCrack  =  hb2n  Equation  3.1 1 

Where  h  [L]  is  the  crack  depth  and  b  [L]  is  the  crack  radius.  The  surface  area  of  the 
crack,  where  diffusion  will  occur,  is  given  by  Equation  3.12: 

SACrack  =  h2bn  Equation  3.12 

where  SA  is  the  crack  surface  area  [L2].  The  volume  of  NAPL  in  a  crack  is  given  by 
Equation  3.13: 

K NAPL  V crack  (1  —  rsw )  Equation  3.13 

Where  Vnapl  [L3]  is  the  volume  of  NAPL  in  the  crack  and  rsw  [-]  is  the  residual  saturation 
of  water.  The  water  residual  saturation  is  the  fraction  of  water  in  the  crack  which  will  not 
be  displaced  by  DNAPL.  The  residual  water  saturation  can  range  from  0  to  1 ,  and  can  be 
determined  for  a  particular  media  based  on  relative  permeability  curves.  In  this  thesis, 
the  residual  water  saturation  is  assumed  to  be  0.1. 

Within  the  model,  modeling  individual  cracks  proved  infeasible;  therefore,  as 
noted  earlier,  flow  in  cracks  was  simulated  as  increased  vertical  conductivity.  In  order 
to  model  the  DNAPL  within  the  cracks  it  is  assumed  that  the  DNAPL  fully  penetrates  the 
cracks  instantaneously.  With  this  assumption  the  total  volume  of  DNAPL  entering  the 
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cracks  can  be  calculated  and  distributed  within  the  matrix.  Equation  3.14  gives  the  total 
volume  of  DNAPL. 

V NAPLtotai  =  Vnaplc  Equation  3.14 

The  DNAPL  will  be  simulated  as  being  evenly  distributed  throughout  the  areas 
containing  cracks  with  saturation  Snapl,  given  below  in  Equation  3.15: 

VcracKsourceSNAPL  =  Equation  3.15 

Where  Vcracksource  [L  ]  is  the  volume  of  the  aquifer  containing  cracks,  Snapl  [-]  is  the 
DNAPL  saturation  (volume  of  DNAPL  per  volume  of  void).  Figure  3.4  is  a  simple 
conceptual  model  of  DNAPL  distribution  at  a  residual  saturation.  The  DNAPL  is  shown 
in  red,  the  soil  grains  in  dark  blue,  and  the  water  in  light  blue. 


DNAPL 


Figure  3-D:  DNAPL  in  Cracks  Approximated  as  Saturation 
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3.4  Model  Implementation 

The  effect  of  cracking  was  examined  using  a  program  called  Ground  Water 
Modeling  System  (GMS).  The  scenarios  discussed  below  are  hypothesized  to 
demonstrate  the  impact  of  cracking  on  down  gradient  contaminant  plume  concentrations. 
The  contaminant  modeled  in  these  simulations  is  TCE. 

3.4.1  Model  Scenarios 

Three  distinct  scenarios  are  considered  and  evaluated  using  GMS.  The  three 
scenarios  are  evaluated  using  the  same  baseline  conditions  for  hydraulic  gradient, 
DNAPL  pool  position,  pool  area,  source  exposure/removal  monitoring  time,  and 
monitoring  well  position.  The  three  scenarios  are  listed  below  and  will  be  referred  to  by 
scenario  number  for  the  remainder  of  this  thesis. 

1 .  Transport  of  the  CAH  into  an  uncracked  clay  matrix  (transport  assumed  to  be 
governed  by  diffusion  only) 

2.  Transport  of  the  dissolved  CAH  into  a  cracked  clay  matrix  (transport  assumed  to 
be  governed  by  advection  and  diffusion) 

3.  Transport  of  the  DNAPL  and  diffusion  of  the  CAH  into  a  cracked  clay  matrix 
(transport  assumed  to  be  governed  by  diffusion  and  advection  of  the  CAH 
coupled  with  dissolution  of  the  DNAPL) 

Figure  3.5  provides  a  conceptual  model  of  each  scenario  beneath  the  source  zone.  Only 
one  crack  is  shown  for  scenarios  2  and  3. 
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(1)  {2)  (3) 


Figure  3-E:  Conceptual  Diagram  of  Contaminant  Transport  Beneath  a  Source  Zone  for  Scenarios  1- 

3 

3.4.1  Ground  Water  Modeling  System 

GMS  was  used  to  evaluate  the  impact  of  cracking  on  down  gradient  plume 
concentrations.  Three  GMS  modeling  packages  were  used  in  this  work,  MODFLOW, 
MODPATH,  and  RT3D.  A  three  dimensional  grid  was  created  with  a  length  of  100m,  a 
width  of  70m,  and  a  total  depth  of  14m.  Each  cell  in  the  grid  was  a  1  meter  cube.  The 
high  permeability  sand  layer  was  6m  thick  atop  a  6m  low  permeability  clay  layer.  The 
clay  layer  overlaid  a  2m  high  permeability  sand  layer.  The  2m  sand  layer  was  needed  to 
avoid  having  a  no-flow  boundary  condition  at  the  bottom  of  the  clay  layer.  Within  the 
low  permeability  layer,  a  192m  cracked  clay  zone  was  emplaced  for  scenarios  2  and  3. 
For  more  detail  on  the  numerical  domain,  see  Appendix  A.  Figure  3.6  is  a  conceptual 
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diagram  depicting  the  media  properties  used  in  scenario  one  (Figure  3.6(a))  and  in 
scenarios  2  and  3  (Figure  3.6(b)).  The  darkened  area  in  Figure  3.6b  represents  the 
cracked  area  beneath  the  source. 

GMS  was  used  to  evaluate  the  impact  of  cracking  on  down  gradient  plume 
concentrations.  Three  GMS  modeling  packages  were  used  in  this  work,  MODFLOW, 
MODPATH,  and  RT3D.  A  three  dimensional  grid  was  created  with  a  length  of  100m,  a 
width  of  70m,  and  a  total  depth  of  14m.  Each  cell  in  the  grid  was  a  1  meter  cube.  The 
high  permeability  sand  layer  was  6m  thick  atop  a  6m  low  permeability  clay  layer.  The 
clay  layer  overlaid  a  2m  high  permeability  sand  layer.  The  2m  sand  layer  was  needed  to 
avoid  having  a  no-flow  boundary  condition  at  the  bottom  of  the  clay  layer.  Within  the 
low  permeability  layer,  a  192m2  cracked  clay  zone  was  emplaced  for  scenarios  2  and  3. 
For  more  detail  on  the  numerical  domain,  see  Appendix  A.  Figure  3.6  is  a  conceptual 
diagram  depicting  the  media  properties  used  in  scenario  one  (Figure  3.6(a))  and  in 
scenarios  2  and  3  (Figure  3.6(b)).  The  darkened  area  in  Figure  3.6b  represents  the 
cracked  area  beneath  the  source. 
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Figure  3-F:  Model  used  for  (a)  Scenario  1  -  uncracked  clay  (b)  Scenarios  2  and  3  -  cracked  clay 
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As  a  first  step  in  modeling  contaminant  transport,  initial  flow  conditions  must  be 


established.  The  aquifer  was  assumed  to  be  unconfined,  and  a  horizontal  hydraulic 
gradient  of  lm/lOOm  was  created  by  setting  the  head  at  the  left  and  right  side  boundaries 
of  the  domain  at  14.5  and  13.5m,  respectively.  A  vertical  hydraulic  gradient  of 
0.5m/14m  was  established  by  setting  the  bottom  heads  at  the  left  and  right  hand 
boundaries  to  14.0m  and  13.0m,  respectively.  Values  of  the  vertical  and  horizontal 
hydraulic  conductivity,  longitudinal  dispersivity,  and  porosity  were  assigned  to  the  sand, 
clay,  and  cracked  clay  as  shown  below  in  Table  3.1. 


Table  3-1:  Parameter  Values  Input  into  MODFLOW 


Sand 

Clay 

Cracked  Clay 

Horizontal  Kh 
(m/d) 

17.28 

4.32xl0"5 

4.32xl0"5 

Vertical  Kv  (m/d) 

1.728 

4.32xl0'6 

4.32x1 0'4 

Longitudinal 
Dispersivity  (m) 

1.0 

l.OxlO'4 

l.OxlO'4 

Porosity  (0) 

0.35 

0.43 

0.487 

Miniter  et  al.  (2011) 


MODFLOW  was  used  to  compute  the  steady  state  flow  heads  for  each  of  the 
three  scenarios  based  on  the  system  parameters  in  Table  3-1.  Three  monitoring  wells 
were  placed  50m  down  gradient  of  the  source  to  monitor  plume  concentrations.  The 
locations  of  the  monitoring  wells  are  shown  below  in  Figure  3.7  as  black  dots.  One  well 
was  placed  in  the  sand  layer,  one  at  the  clay/sand  interface,  and  one  in  the  clay.  A  typical 
MODFLOW  output  showing  hydraulic  heads  throughout  the  model  domain  is  shown 
below  in  Figure  3.7. 
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Figure  3-G:  MODFLOW  Output  Along  a  Longitudinal  Cross  Section 

MODPATH  can  then  be  used  to  evaluate  the  flow  of  particles  through  the  system  to  help 
visualize  flow  conditions  and  identify  stagnation  points. 

The  MODFLOW  solution  is  then  read  by  RT3D  which  then  simulates  the 
advection,  dispersion,  sorption,  and  diffusion  of  the  contaminant  through  the  domain. 
Two  periods  were  used  to  simulate  a  DNAPL  release  into  the  subsurface.  In  the  first  ten 
year  period,  a  1 92m2  DNAPL  pool  source  was  simulated  by  holding  the  concentration  of 
cells  at  a  constant  value  for  10  years.  After  this  ten  year  period,  it  was  assumed  the 
source  is  remediated,  and  the  constant  concentration  source  was  removed.  The 
contaminants  were  then  transported  out  of  the  low  permeability  layer  based  on 
concentration  gradients.  The  model  was  then  run  for  an  additional  40  years  to  examine 
the  effect  of  this  back  diffusion  on  down  gradient  plume  concentrations  as  well  as  the 
DNAPL  dissolution  into  the  aquifer. 

In  RT3D  a  user  defined  transport  package  written  by  Dr.  Junqi  Huang  (Personal 
communication,  2011)  was  used.  The  transport  package  requires  the  user  to  define 
various  parameter  values  for  the  contaminant  and  media.  These  parameters  are  the  bulk 
density  of  the  media,  the  sorption  coefficient  of  the  media,  the  mass  transfer  coefficient 
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from  the  dissolved  phase  to  the  soil,  solubility  of  the  DNAPL,  bulk  density  of  the 
DNAPL,  the  mass  transfer  coefficient  describing  the  dissolution  of  the  DNAPL  to  the 
surrounding  water,  and  the  diffusion  coefficient  in  the  media.  The  values  used  for  the 
different  media  are  shown  below  in  Table  3.2  and  the  contaminant  properties  are  shown 
below  in  Table  3.3. 


Table  3-2:  RT3D  Media  Values 


Media 

Sand 

Clay 

Cracked  Clay 

Bulk  Density 
(pb)  (kg-m"3) 

1722 

1499 

1349 

Sorption 
Constant  (Kd) 
(mg-L1) 

65.4 

1908 

1908 

Diffusion 
Coefficient 
(Deff)  (cm2-s  -1) 

0.0 

8.64  x  1  O'6 

8.64  x  1  O’6 

(Miniter,  20] 

H) 

Table  3-3:  RT3D  Chemical  Values 


0>i 


Media 

TCE 

Bulk  Density  (pb)  (kg-m'3) 

1460 

NAPL  Solubility  (Cs)  (mg-L1) 

110(1) 

First  Order  Sorption  Rate 
Constant  (a)  (d1) 

0.0 

First  Order  NAPL  Dissolution 
Rate  Constant  (P)  (d'1) 

0.00445(2) 

2006) 


Sorption  is  not  examined  in  this  work.work;  therefore  the  first  order  sorption  rate 
constant  equal  used  is  zero,  effectively  removing  all  sorption  processes.  Using  Equations 
3.14  and  3.15  the  NAPL  saturation  is  determined  for  scenario  3.  For  these  scenarios  it  is 
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assumed  the  average  crack  aperture  is  150  pm,  the  cracks  are  0.03m  apart,  and  extend  4m 
in  depth.  Given  these  values  it  is  assumed  there  are  214,401  cracks  in  the  192  m2  source 
zone.  The  DNAPL  saturation  is  set  at  6.62x1  O'5  for  scenario  3  at  two  times,  0  and  ten 
years.  The  saturation  cannot  be  held  constant  throughout  the  first  ten  years,  therefore  it  is 
reset  at  ten  years  when  the  source  has  been  removed  to  more  realistically  simulate 
DNAPL  recharge  into  the  cracks. 

3.5  Results  Analysis 

Two  methods  were  used  to  quantify  the  effect  of  cracks  on  DNAPL  contaminant 
fate  and  transport:  a  mass  analysis  and  examination  of  concentration  versus  time 
breakthroughtime  breakthrough  curves  at  down  gradient  locations.  The  first  method,  a 
mass  analysis,  quantified  the  total  mass  stored  in  the  low  permeability  layer.  Mass  could 
be  stored  in  three  ways:  in  the  aqueous  phase,  sorbed  to  the  soil  solids  (though  in  this 
study,  sorption  is  considered  negligible),  and  in  Scenario  3,  in  the  NAPL  phase.  The 
mass  in  the  three  scenarios  was  compared  at  three  points  in  time:  at  the  start  of  the 
simulation,  immediately  after  the  source  is  removed,  and  40  years  after  the  source  has 
been  removed. 

The  second  method,  used  by  Miniter  et  al.  (2011),  examined  predicted 
concentrations  at  the  observation  wells.  Typically,  it  will  take  a  relatively  short  time  for 
dissolved  contaminant  from  the  source  to  reach  the  observation  well.  After  the  source  is 
removed,  the  concentration  at  a  given  observation  well  will  decline.  It  has  been  observed 
that,  with  DNAPLs,  the  concentration  rapidly  declines  when  the  source  is  removed, 
followed  by  the  persistence  of  low  levels  for  extended  periods  of  time.  This  slow  decline 
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at  long  periods  of  time  is  termed  “tailing”  (Parker  et  al.,  2008).  This  concentration 
history  can  then  be  used  to  determine  the  time  required  for  aqueous  concentrations  to 
reach  the  regulatory  MCL  (5  ug/L  for  TCE). 

It  is  also  possible  to  describe  the  breakthrough  curve  by  its  first  moment.  The 
first  moment  measures  the  center  of  mass  of  a  distribution.  The  higher  the  value  of  the 
first  moment  of  a  breakthrough  curve  at  a  given  monitoring  well,  the  longer 
contamination  is  persisting  at  that  monitoring  well.  A  high  value  for  the  first  moment 
may  be  indicative  of  tailing.  It  is  hypothesized  that  both  the  first  moment  and  the  time  to 
attain  MCLs  will  increase  significantly  from  Scenario  1  to  3. 


In  this  study,  the  first  moment  of  breakthrough  curves  was  determined  for 
Scenarios  1-3. 

In  order  to  calculate  the  first  moment,  the  area  under  the  breakthrough  curve  was 
calculated  using  Equation  3.16. 

area  =  Z|™oX_1  p^tl+1^+c(-tl')j  (ti+1  —  tj)  Equation  3.16 

where  c(t)  [mg-L1]  is  the  concentration  output  from  GMS  at  time  t,  [T], 

The  discrete  residence  time  density  function  (f(t;))  is  then  determined  for  each  discrete 
concentration  value  using  Equation  3.17. 

f(tj)  =  Equation  3.17 


The  first  moment  is  then  calculated  using  Equation  3.18. 

iRTD  =  SILT'1  p^2]  [f(ti)+fi+l)]  (t1+i  -  ti)  Equation  3.18 


where  tRTD  [T]  is  the  mean  residence  time  or  first  moment  (Clark,  2009). 
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3.5.1  Sensitivity  Analysis 

A  sensitivity  analysis  was  used  to  determine  those  parameters  which,  when 
varied,  have  the  most  impact  on  down  gradient  plume  concentrations.  Scenario  3  was  be 
used  in  the  simulations  to  quantify  the  impact  of  different  model  parameters.  The 
sensitivity  analysis  focused  on  the  impact  of  varying  the  first  order  NAPL  dissolution  rate 
constant  (P)  and  the  vertical  hydraulic  conductivity  (Kv)  below  the  source.  P  was  varied 
by  two  orders  of  magnitude  in  order  to  determine  the  effect  of  DNAPL  dissolution.  As  P 
increases  or  decreased,  the  concentration  at  the  source  should  also  increase  or  decrease  in 
accordance  with  Equations  3.7  and  3.8.  The  vertical  conductivity,  Kv,  was  also  varied  by 
two  orders  of  magnitude  to  determine  the  impact  of  increased  and  decreased  vertical  flow 
at  the  source  zone.  Table  3.4  below  lists  baseline  values  for  p  and  Kv  as  well  as  the 
values  used  in  this  sensitivity  analysis. 


Table  3-4:  Sensitivity  Analysis  Values 


Lower  Sensitivity 
Analysis 

Vaslue  Value 

Baseline 

Value 

Upper  Sensitivity 
Analysis  Value 

3  (d1) 

4.45E-5 

0.00445 

0.445 

Kv(m-d') 

4.32E-6 

0.000432 

0.0432 

The  effect  of  changing  p  and  Kv  will  be  quantified  by:  (1)  time  to  reach  MCL  at 
down  gradient  wells,  (2)  down  gradient  well  breakthrough  curve  first  moment,  and  (3) 
comparison  of  mass  of  contaminant  stored  in  the  low  permeability  layer  immediately 
following  source  removal  (10  years)  and  at  the  conclusion  of  the  simulation  (50  years). 
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4.  Results  and  Analysis 


4.1  Overview 

This  chapter  presents  the  results  of  the  simulations  described  in  Chapter  3.  These 
results  include  breakthrough  curves  for  Scenarios  1-3,  mass  storage  comparisons,  and  the 
sensitivity  analysis.  The  first  moment  analysis  discussed  in  Chapter  3  was  not  used  due 
to  the  absence  of  significant  differences  in  the  breakthrough  curves. 

4.2  Simulation  Results 

The  figures  and  data  presented  in  this  section  are  for  Scenarios  1-3.  The  MCL  for 
TCE,  0.005  mg-L'1,  is  shown  in  all  plots  for  comparison. 

4.2.1  Breakthrough  Curves 

Three  breakthrough  curves  are  shown  below  in  Figure  4.1.  These  breakthrough 
curves  plot  concentration  versus  time  for  the  total  simulation  time,  50  years  (18250  days). 
Figure  4.  la  is  a  plot  of  the  concentrations  at  the  observation  point  placed  within  the  sand, 
Figure  4.  lb  is  a  plot  of  the  concentrations  observed  at  the  observation  point  placed  at  the 
sand/clay  interface,  and  Figure  4.1c  is  a  plot  of  the  concentrations  observed  at  the 
observation  point  placed  within  the  clay  layer.  All  wells  are  50m  down  gradient  from  the 
source  zone. 
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Scenario  1  - 
Uncracked  Clay 

Scenario  2  -  Cracked 
Clay 

Scenario  3  -  Cracked 
Clay  with  DNAPL 


- MCL 


(b) 


60 


Scenario  1  -  Uncracked 
Clay 

Scenario  2  -  Cracked 
Clay 

Scenario  3  Cracked  Clay 
with  DNAPL 

- MCL 


(C) 

Figure  4-A:  Breakthrough  Curves  (a)  within  sand  layer,  (b)  at  sand/clay  interface,  and  (c)  within  clay 

layer 


Qualitatively,  the  breakthrough  curves  for  Scenarios  1-3  all  display  the  expected 
behavior.  At  the  observation  points  in  high  flow  sand  zones  (4. la  and  4. lb),  the 
concentration  rises  quickly  when  the  source  is  present,  and  when  the  source  is  removed  a 
rapid  decrease  is  observed  with  back  diffusion  then  causing  the  tailing.  In  the  low  flow 
clay  zone,  Figure  4.1c,  the  increase  was  much  slower  due  to  low  flow  rates  in  the  low 
permeability  layer  and  a  relative  increase  in  diffusive  rather  than  advective  transport  into 
the  low  permeability  layer. 


Contrary  to  results  presented  in  Miniter  (2011),  the  concentration  in  the  sand  layer 
does  not  remain  above  the  MCL  for  extended  periods  of  time  in  Scenario  2.  However, 
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consistent  with  the  findings  of  Miniter  (2011),  at  the  sand  clay  interface  and  within  the 
clay  layer,  the  concentrations  remain  above  the  MCL  for  over  forty  years.  Furthermore, 
the  tailing  seen  at  the  sand-clay  interface  is  similar  to  the  tailing  from  back  diffusion  of 
TCE  in  a  clayey  silt  aquitard  observed  by  Chapman  and  Parker  (2005).  To  further 
illustrate  differences  between  the  scenarios,  curves  are  presented  in  Figure  4.2  over  the 
time  period  beginning  ten  years  after  source  removal  (7300  days)  and  extending  to  fifty 
years. 


Scenario  1  - 
Uncracked  Clay 

Scenario  2  -  Cracked 
Clay 

Scenario  3  -  Cracked 
Clay  with  DNAPL 


(a) 
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Scenario  1  - 
Uncracked  Clay 

Scenario  2  -  Cracked 
Clay 

Scenario  3  -  Cracked 
Clay  with  DNAPL 


- MCL 


(b) 


Scenario  1  - 
Uncracked  Clay 

Scenario  2  -  Cracked 
Clay 

Scenario  3  Cracked 
Clay  with  DNAPL 


- MCL 


(C) 

Figure  4-B:  Breakthrough  Curves  (a)  within  sand  layer,  (b)  at  sand/clay  interface,  and  (c)  within  clay 

layer 
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As  Figure  4.2a  shows,  the  concentrations  of  TCE  within  the  sand  layer  and  clay 
for  all  the  scenarios  appear  to  be  the  same  for  the  entire  duration  of  the  simulation.  These 
concentrations  remain  well  below  the  MCL,  indicating  cracking  may  have  little  to  no 
impact  on  down  gradient  concentrations  1  m  above  the  low  permeability  layer.  However, 
the  concentrations  observed  at  the  sand-clay  interface,  shown  in  Figure  4.2b,  do  indicate 
a  significant  difference  in  down  gradient  concentrations  associated  with  no  cracking, 
cracking,  and  DNAPL  storage  and  transport  within  cracks. 

While  the  concentrations  at  the  sand-clay  interface  remain  within  an  order  of 
magnitude  of  each  other,  the  time  at  which  the  concentration  decreases  below  the  MCL  is 
significantly  different.  In  Scenario  1  the  concentration  drops  below  the  MCL  at  40.3 
years,  in  Scenario  2  the  corresponding  time  is  44.6  years,  and  in  Scenario  3  the 
concentration  remains  above  the  MCL  for  the  extent  of  the  simulation.  Further 
simulations  concluded  Scenario  3  dropped  below  the  MCL  at  52.6  years.  Since  the 
concentration  breakthrough  curves  within  the  clay  layer,  shown  in  Figure  4.2c,  are  all  the 
same,  the  differences  observed  at  the  sand/clay  interface  are  most  likely  not  due  to  back 
diffusion  in  the  vicinity  of  the  observation  point,  but  rather  due  to  up  gradient 
contamination.  This  behavior  can  be  explained  by  the  liquid  phase  TCE  stored  within  the 
cracks. 

4.2.2  Mass  Analysis 

The  goal  of  the  mass  analysis  is  to  quantify  the  amount  of  TCE  stored  within  the 
low  permeability  layer  at  different  points  in  time;  at  the  start  of  the  simulation  (0  years), 
immediately  following  source  removal  (10  years),  and  at  the  end  of  the  simulation  (50 
years).  The  mass  of  stored  TCE  is  the  long  term  source  for  back  diffusion,  and  therefore, 
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an  important  parameter  to  measure.  The  total  mass  is  calculated  for  the  entire  100  m  x  70 
m  x  6  m  clay  layer,  although  the  mass  distribution  is  not  determined.  It  is  believed, 
though,  that  much  of  the  mass  is  in  the  upper  portion  of  the  clay  layer.  It  is  important  to 
note  in  order  to  simulate  the  effect  of  cracks  in  Scenario  3,  the  DNAPL  saturation  was  set 
at  the  beginning  of  the  simulation  and  again  at  the  ten  year  point,  as  discussed  in  Chapter 
3.  Unfortunately  this  means  between  0  and  10  years  the  DNAPL  saturation  is  not  held 
constant.  Table  4.1  contains  calculated  TCE  masses  in  the  clay  layer  for  Scenarios  1-3. 


Table  4-1:  Mass  Analysis  for  Scenarios  1-3 


Scenario  1  Scenario  2  Scenario  3 


0  yrs 

10  yrs 

50  yrs 

0  yrs 

10  yrs 

50  yrs 

0  yrs 

10  yrs 

50  yrs 

Dissolved  TCE  (kg) 

0 

8.19 

6.23 

0 

8.29 

6.36 

0 

8.38 

6.89 

Sorbed  TCE  (kg) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

DNAPL  TCE  (kg) 

0 

0 

0 

0 

0 

0 

36.1 

36.1 

35.7 

Total  Mass  TCE  (kg) 

0 

8.19 

6.23 

0 

8.29 

6.36 

36.1 

44.4 

42.6 

As  expected,  the  dissolved  TCE  mass  at  both  10  and  50  years  increased  from 
Scenarios  1  to  3.  While  the  increase  in  mass  is  not  significant,  it  may  explain  the  tailing 
seen  at  the  sand-clay  interface  in  Figure  4.2b.  In  both  Scenarios  1  and  2  the  dissolved 
mass  decreased  by  approximately  1.9  kg  between  years  10  and  50,  while  in  Scenario  3 
the  decrease  was  only  1 .5  kg.  The  higher  aqueous  concentrations  observed  in  Scenario  3 
are  most  likely  due  to  the  dissolution  of  the  DNAPL.  Further  examination  is  required  to 
determine  if  dissolution  of  TCE  will  sustain  aqueous  concentrations  in  the  aquifer  above 
the  MCL. 


4.4.3  Sensitivity  Analysis 

The  goal  of  the  sensitivity  analysis  was  to  determine  which  input  values 
significantly  impacted  breakthrough  tailing  at  the  observation  points  and  the  total  mass 
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storage  in  the  low  permeability  layer.  Scenario  3  was  used  as  the  baseline  scenario  for 
this  analysis  in  order  to  capture  all  important  processes.  The  values  changed  include  the 
first  order  DNAPL  dissolution  rate  constant  (P),  and  the  vertical  hydraulic  conductivity 
(Kv)  of  the  cracked  clay  directly  below  the  source.  Table  4.2  indicates  both  the  baseline 
and  sensitivity  analysis  values. 


Table  4-2:  Sensitivity  Analysis  Values 


Lower  Sensitivity 
Analysis  Vaslue 

Baseline 

Value 

Upper  Sensitivity 
Analysis  Value 

P  (d ') 

4.45E-5 

0.00445 

0.445 

Kv(m-d_1) 

4.32E-6 

0.000432 

0.0432 

With  the  changes  in  p  and  Kv  no  discernible  differences  were  observed  in  the  clay  and 
sand  breakthrough  curves;  however,  noticeable  changes  occurred  at  the  sand-clay 
interface.  The  breakthrough  curve  at  the  sand-clay  interface  for  varying  p  and  Kv  is 
presented  below  in  Figure  4.3(a)  and  (b).  The  plot  includes  the  curves  generated  using 
the  baseline  values  as  well  as  using  the  sensitivity  analysis  values.  A  mass  analysis  was 
used  in  conjunction  with  breakthrough  curves  to  quantify  the  effect  of  varying  P  and  Kv. 
The  results  of  the  mass  analysis  are  shown  below  in  Table  4.3.  The  results  of  Scenario  3 
with  baseline  values  are  also  shown  in  Table  4.3  for  comparison. 


66 


Baseline  -  Cracked 
Clay  with  DNAPL 

Vertical 

Conductivity=0 .0432 
(m/d) 

—Vertical 

Conductivity=4 . 3  2E-6 
(m/d) 


(a) 


(b) 


Figure  4-C:  Breakthrough  Curve  at  Sand-Clay  Interface  for  Sensitivity  Analysis:  (a)  p  and  (b)  Kv 
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Table  4-3:  Mass  Analysis  for  Sensitivity  Analysis 


Scenario  3  -  p=0. 0000455  d"1  Scenario  3  Scenario  3  -  p=0.455  d" 


0  yrs 

10  yrs 

50  yrs 

0  yrs 

10  yrs 

50  yrs 

0  yrs 

10  yrs 

50  yrs 

Dissolved  TCE  (kg) 

0 

8.3 

6.5 

0 

8.38 

6.89 

0 

15.6 

29.4 

Sorbed  TCE  (kg) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

DNAPL  TCE  (kg) 

36.1 

36.1 

36.1 

36.1 

36.1 

35.7 

36.1 

28.7 

0.1 

Total  Mass  TCE  (kg) 

36.1 

44.3 

42.6 

36.1 

44.4 

42.6 

36.1 

44.3 

29.5 

Scenario  3 

1  -  Kv=4.32E-6  m-d"1 

Scenario  3 

Scenario  ! 

3  -  Kv=0.0432  m-d-1 

0  yrs 

10  yrs 

50  yrs 

0  yrs 

10  yrs 

50  yrs 

0  yrs 

10  yrs 

50  yrs 

Dissolved  TCE  (kg) 

0 

8.22 

6.73 

0 

8.38 

6.89 

0 

8.61 

7.13 

Sorbed  TCE  (kg) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

DNAPL  TCE  (kg) 

36.1 

36.1 

35.8 

36.1 

36.1 

35.7 

36.1 

36.1 

35.8 

Total  Mass  TCE  (kg) 

36.1 

44.3 

42.5 

36.1 

44.4 

42.6 

36.1 

44.7 

42.9 

When  p  was  increased,  a  large  increase  in  tailing  was  noticed.  As  P  is  the 
parameter  that  represents  the  effects  of  dissolution,  the  increase  in  tailing  can  be 
attributed  to  the  dissolution  of  the  DNAPL.  Since  the  rate  of  dissolution  is  larger,  the  rate 
of  mass  transport  of  CAH  into  the  aqueous  phase  and  into  the  high  permeability  layer  is 
larger.  Note  that  due  to  model  limitations,  the  extent  of  tailing  may  be  even  more 
pronounced  in  a  real  scenario,  as  the  model  used  in  this  study  does  not  allow  for 
replenishment  of  the  DNAPL  in  the  cracks  by  the  DNAPL  pool  sitting  atop  the  clay 
layer.  As  Figure  4.3a  shows,  it  appears  that  the  concentration  is  still  increasing  at  the  end 
of  the  simulation.  This  trend  is  also  observed  in  the  mass  analysis,  which  shows  the 
aqueous  concentration  at  the  source  increasing  significantly  after  the  source  is  removed, 
due  to  faster  DNAPL  dissolution  at  the  source.  This  indicates  DNAPL  within  the  cracks 
at  the  source  can  contribute  significantly  to  down  gradient  plume  concentrations  for  an 
extended  period  of  time.  When  p  was  decreased,  no  change  was  detected  from  the 
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baseline  scenario.  This  indicates  the  value  used  in  Scenario  3  may  be  artificially  low.  We 
should  note  that  the  p  values  employed  in  these  simulations  were  based  on  pure  TCE,  not 
TCE  waste  found  at  contaminated  sites.  A  quantification  of  p  for  TCE  waste  is  needed  to 
realistically  simulate  the  behavior  of  DNAPL  waste  at  sites  where  the  DNAPL  waste  has 
penetrated  cracks. 

In  this  research,  cracks  are  approximated  by  an  equivalent  hydraulic  conductivity; 
as  cracking  increases  so  does  the  vertical  hydraulic  conductivity  (Kv).  The  vertical 
conductivity  was  increased  to  account  for  cracks  based  on  the  general  observation  that 
cracking  increases  the  overall  vertical  hydraulic  conductivity  by  two  to  three  orders  of 
magnitude.  A  correlation  between  vertical  hydraulic  conductivity  and  crack 
characteristics  such  as  crack  aperture,  spacing,  and  depth  would  be  useful  in  determining 
reasonable  estimates  for  vertical  hydraulic  conductivity  in  future  simulations.  In  this 
sensitivity  analysis,  when  the  vertical  conductivity  increased,  tailing  also  increased. 

After  the  initial  ten  years,  the  mass  within  the  clay  layer  in  the  Scenario  3  with  the 
increased  vertical  hydraulic  conductivity  is  similar  to  that  in  the  baseline  Scenario  3; 
however,  40  years  after  source  removal,  the  mass  in  the  former  scenario  is  higher  than 
that  in  the  latter.  The  results  of  this  sensitivity  analysis  indicate  more  cracking  can  lead 
to  higher  down  gradient  concentrations  for  extended  periods  of  time. 
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5  Conclusions  and  Recommendations 


5.1  Conclusions 

This  research  modeled  the  effect  of  cracking  in  low  permeability  layers  on  the 
storage  and  transport  of  TCE.  The  effect  of  cracking  on  transport  was  measured  by 
comparing  breakthrough  curves  at  down  gradient  monitoring  points  as  well  as  by 
calculating  the  quantity  of  TCE  stored  in  the  low  permeability  clay  layer. 

The  model  simulated  cracks  through  an  increase  in  the  vertical  conductivity  at  the 
source.  This  change  was  hypothesized  to  simulate  enhanced  diffusion  of  TCE  into  the 
low  permeability  layer.  In  the  model  the  source  was  present  for  ten  years  allowing  for 
contaminant  transport  into  and  storage  within  the  low  permeability  layer.  The  source  was 
then  removed  and  down  gradient  concentrations  were  simulated  at  three  observation 
points  for  an  additional  forty  years  to  determine  the  impact  of  cracking.  Three  scenarios 
were  studied,  (1)  transport  into  uncracked  clay,  (2)  transport  of  the  dissolved  TCE  into 
cracked  clay,  and  (3)  transport  of  DNAPL  phase  TCE  into  cracked  clay.  Down  gradient 
concentrations  were  sustained  by  back  diffusion  of  TCE  from  the  low  permeability  layer. 
Based  on  the  results  of  the  mass  analysis  and  breakthrough  curves  it  was  determined  that: 

(1)  Cracking  (as  approximated  with  an  increase  in  vertical  conductivity  and  a  DNAPL 
saturation)  will  cause  an  increase  in  TCE  transport  into  the  low  permeability  layer 

(2)  Enhanced  transport  of  TCE  into  the  source  zone  will  sustain  down  gradient 
concentrations  above  the  MCL  at  the  sand-clay  interface 
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(3)  Down  gradient  concentrations  are  sustained  due  to  back  diffusion  from  the  source 


zone 

(4)  DNAPL  phase  TCE  within  cracks  can  significantly  contribute  to  down  gradient 
concentrations;  however,  this  contribution  is  dependent  upon  the  rate  of  DNAPL 
dissolution 

(5)  Remediation  goals  may  be  impossible  to  meet  within  a  prescribed  time  frame  if 
source  remediation  strategies  are  used  which  do  not  account  for  the  increased  back 
diffusion  out  of  cracked  low  permeability  layers  at  contaminated  sites. 

5.2  Recommendations  for  Future  Research 

The  limitations  of  the  model  used  in  this  research  inherently  provide  future 
opportunities  for  research.  This  work  approximates  cracking  in  the  source  zone  through 
an  increase  of  the  vertical  hydraulic  conductivity.  While  this  assumption  may  permit 
enhanced  transport,  the  use  of  an  adaptable  grid  model  would  allow  for  the  direct 
modeling  of  cracks.  If  an  adaptable  grid  model  is  used,  it  might  be  possible  to  model 
pressure  dependent  NAPL  transport  into  the  cracks  for  the  duration  of  source  exposure. 
Further,  direct  crack  modeling  would  only  simulate  increased  flow  through  cracks  as 
opposed  to  the  whole  matrix,  meaning  diffusion  out  of  the  uncracked  clay  may  in  fact  be 
a  much  slower  process.  The  results  of  Ayral  et  al.  (2011)  indicate  that  TCE  DNAPL 
waste  can  cause  cracking  in  clay;  this  means  if  DNAPL  is  stored  within  existing  cracks, 
the  crack  properties  may  change  allowing  for  further  diffusion  into  the  matrix. 

In  order  to  evaluate  the  ability  of  this  model  to  predict  field  scale  conditions,  the  model 
should  be  applied  to  a  well  characterized  site  that  is  known  to  display  back  diffusion  after 
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source  zone  remediation.  While  cracking  may  be  unknown  at  the  site,  a  variety  of 
hydraulic  conductivities  can  be  tested  to  determine  if  cracking  can  be  an  explanation  of 
long  term  plume  persistence  above  target  MCLs. 

This  research  used  dissolution  and  diffusion  values  for  pure  TCE.  TCE  waste  can  have 
very  different  properties  due  to  the  presence  of  other  chemicals;  therefore  further  research 
could  determine  and  incorporate  realistic  chemical  properties  of  TCE  waste  into 
simulations.  This  model  also  does  not  account  for  chemical  or  biological  contaminant 
degradation.  TCE  degradation  could  significantly  impact  the  down  gradient 
concentration;  however,  if  degradation  is  modeled  it  would  be  important  to  also  model 
TCE  daughter  product  formation. 
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Appendix  A 

Open  GMS 

Start  Up 

Save  your  file  -  Be  sure  to  save  after  every  few  steps  as  GMS  can  crash  frequently. 

The  units  are  not  assumed  to  be  SI  units.  The  units  used  in  the  Sievers  thesis  are  SI  units. 
Edit  =>  Units... 

The  required  packages  for  running  simulations  are  MODFLOW  and  MT3D,  RT3D, 
SEAM3D. 

Edit  =>  Model  Interfaces  =>  Check  MODFLOW  and  MT3D,  RT3D,  SEAM3D 

boxes 

Create  Grid 

On  Bottom  of  screen  click  3D  Grid  symbol  (Green  3D  cube)  shown  in  Figure  A.  1 


Figure  A.  1:  3D  grid  symbol 


Once  the  3D  Grid  is  selected  a  Grid  dropdown  menu  will  appear  in  the  upper  toolbar  next 
to  Display 

Grid  =>  Create  Grid 

Input  dimensions  of  grid.  Make  origin  at  (0,0,0).  Sievers  grid  is  x  =  100m  y=70m 
z=12m.  Be  sure  to  put  number  of  cells  the  same  as  dimensions  to  ensure  lm3  grid  blocks. 
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Input  Matrix  Characteristics 


Edit  =>  Materials  —  See  table  3.1  for  relevant  characteristics 

Create  MODFLOW  Simulation 

MODFLOW  =>  New  Simulation 

The  MODFLOW  Global/Basic  Package  will  pop  up 


•  MODFLOW  Version  -  Use  MODFLOW  2000  version  (2005  version  not 
compatible  with  Huang  model) 

•  Run  Options  -  Select  Forward  Run 

•  Model  Type  -  Select  Steady  State 

•  Packages  =>  Flow  Package  -  Select  Layer  Property  Flow  (LPF) 

•  Packages  =>  Solver  -  Pre-Cond.  Conj.  Grad  (PCG2) 

•  IBOUND  =>  In  all  layers  set  farthest  left  and  right  column  to  -1,  all  other  cells  to 
1  shown  in  Figure  A.2 


r 


E  □  .  •  ai  *  % 


•  i  «*a am*  • 


L, 


Figure  A.  2:  IBOUND  input 
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•  Starting  Heads  =>  The  only  cells  which  will  remain  constant  are  the  ones  with  a  - 
1  value  in  the  IBOUND  array.  Sievers  thesis  set  1  m  head  difference,  farthest  left 
column  to  20.5  farthest  right  to  19.5m  shown  in  Figure  A.3. 


Figure  A.  3:  Initial  head  conditions 

•  Top  Elevation/Bottom  Elevation  =>  Sanity  check,  should  be  lm  difference 
between  each  layer,  as  well  as  a  lm  difference  for  top  and  bottom  for  same  layer 

MODFLOW  =>  LPF  Package 

•  Layer  Property  Entry  Method  -  Select  use  material  IDs  (seems  to  be  easiest  way) 

•  Layer  type  -  Ensure  for  Layer  1  confined  is  selected 

•  Vertical  hydraulic  conductivity  -  Specify  Kv  for  all  layers 

•  Material  IDs  =>  For  grid  input  numbers  corresponding  to  material  type  (Layers  1- 
6  are  sand,  Layers  7-12  are  a  combination  of  cracked  and  uncracked  clay) 

•  Interblock  transmissivity  -  Harmonic  mean 

•  Cell  wetting  parameters  -  Do  not  select  Allow  wetting  of  cells 
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MODFLOW  =>  Run  MODFLOW 


Create  Observation  Wells 

Right  click  on  Project  in  Project  Explorer  Window  =>  New  =>  Conceptual  Model  shown 
in  Figure  A.4 


Project  Explorer 
H  Project 

i..0  @  3D  Grid  Data 
!  S-EIH  grid 

i-a  £CS3  (MODFLOW) 
ta  SCS2  (MODFLOW) 


Figure  A.  4:  Image  of  project  explorer 


Conceptual  Model  Properties  box  will  pop  up  -  Select  Transport  -  Select  RT3D  -  Select 

User  Defined  Reaction  -  Define  aqueous  species  name  -  Click  OK 

Map  Data  will  now  be  in  Project  Explorer  Window  -  Right  click  New  Model  -  select 

New  Coverage 

Coverage  Setup  will  pop  up  -  Select  defined  species  under  Observation  Points  -  Click 
OK  shown  in  Figure  A.5 
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iew  coverage 


Coverage  Setup 


Coverage  name: 


Coverage  type: 


I 


Sounces/Sinks/BCs 


Source/Sink/BC  Type 

A- 

r 

All 

r 

Layer  range 

r 

Wells 

— 

r 

Wells  (MNW) 

r 

Refine  points 

r 

Specified  Head  (CHD} 

r 

Specified  Head  (1  BOUND) 

r 

Specified  How 

•w 

i?  i 

m  ~l !  H 

Default  layer  range:  1  to 


Areal  Properties 


Property  ^ 

r 

lf| 

r 

Color 

r 

Layer  range 

r 

Recharge  rate 

r 

Recharge  cone. 

r 

Horizontal  K 

r 

Vertical  K 

r 

Horizontal  anis. 

* 

rrr  |  t 

Default  elevation: 


Observation  Points 


0.0 


Use  to  define  model  boundary  (active  area] 


3D  grid  layer  option  for  obs.  pts.: 

By  z  location 

T 

MO  DAE  M  models: 

NONE 

T 

Help... 


OK 


Cancel 


Figure  A.  5:  Coverage  set  up 

Right  click  New  Coverage  -  Attribute  Table  -  Enter  coordinates  for  observation  wells  - 
select  obs.  pt  for  Type  shown  in  Figure  A.6 
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3  Properties 


I  XS  I 


Points 

Show: 

All 

BC  type: 

ANY/NONE  w 

Show  point  coordinates 


ID 

Name 

X 

Y 

z 

Type 

Obs. 

COA 

Obs.  COA 
interval 

Obs.  COA 
conf{%) 

Obs.  ( 
std.  i 

All 

w 

1 

downclay 

85.0 

350 

3.0 

obs.  pt 

■W 

1.0 

95 

0.51021 

2 

downclaysand 

85.0 

35.0 

9.0 

obs.  pt 

w 

w 

1.0 

95 

0.51921 

3 

downsand 

85.0 

35.0 

10.0 

obs.  pt 

1.0 

95 

0.51021 

Help... 


Add  Point 


Delete  Point 


OK 


Cancel 


Figure  A.  6:  Observation  point  properties 

Create  RT3D  Simulation 

MT3D  =>  New  Simulation 

The  MT3D  Basic  Transport  Package  will  pop  up  shown  in  Figure  A.7 

•  Model  -  Select  RT3D 

•  Stress  Periods  -  Can  be  changed  by  user  for  different  simulations 

•  Output  Control  -  Ensure  print  or  save  at  specified  interval  is  chosen  for 
breakthrough  curve  data,  always  select  Save  binary  concentration  file  and  Save 
mass  balance  file 

•  Packages  -  Select  Advection  Package,  Dispersion  Package,  Source/sink  mixing 
package,  Chemical  reaction  package  (select  User-defined  Reaction),  GCG  solver 
package 

•  Define  Species  -  Define  three  species,  COA,  NAPLsoil,  NAPLsaturation.  Order 
does  matter.  COA  is  the  only  mobile  species. 


78 


•  Layer  Data  -  Select  use  materials  for  porosity  and  long.  Dispersivity  and  HTOP 
equals  top  of  layer  1 

•  Starting  Concentration  -  Can  be  defined  for  all  species  based  on  goal  of 
simulation 


Figure  A.  7:  RT3D  Basic  Transport  Package  window 
MT3D  =>  Advection  Package  -  Select  Third  order  TVD  scheme  (ULTIMATE) 

MT3D  =>  Dispersion  Package  —  Input  appropriate  parameters 
MT3D  =>  Source/Sink  Mixing  Package  —  This  is  where  the  simulated  source  is  input 
-  need  one  data  point  for  each  grid  block  (192  for  192m3  source  zone)  shown  below 
in  Figure  A. 8. 

•  Layer  -  Should  be  the  layer  immediately  above  the  low  permeability  layer 
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•  Type  -  specconc 

•  COA  -  This  column  will  be  whatever  you  define  as  your  mobile  species 


§  Source/Sink  Mixing  Package 


S3 


(o)  Use  river  package  Use  stream  package 

Max  number  of  all  point  sounces/sinks  in  flow  model 
MXSS:  16B0 


.Areally  distributed  sounces/sinks 


Recharge... 


NAPLsoil 

NAPLsaturation 


Evapotrans.. 


Initialize  point  sounces/sinks  from  MOD  FLOW 


Con.  Head 


Well 


River/Stream 


Gen  Head 


Reset 


Help.. 


Point  sounces/sinks 


Stress  period:  1 

Start  Time:  0.0 
End  Time:  3650.0 


I  I  Use  previous 


Add 


now 
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type 
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A- 
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6 
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6 
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Figure  A.  8:  Source/Sink  Mixing  Package  initial  conditions 

MT3D  =>  Chemical  Reaction  Package  =>  Define  Parameters  (Order  matters,  also,  cannot 

be  spatially  varied)  shown  in  Figure  A.9 
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•  Modop  -  Always  select  6 

•  Rhobulk  -  Bulk  density  of  media  (mg/L)  -  Note,  only  one  bulk  density 
can  be  set  since  it  cannot  be  spatially  varied,  so  for  entire  array  use  clay 
bulk  density. 

•  Rho  NAPL  -  Bulk  density  of  NAPL  (mg/L) 

•  Alpha  -  First  order  mass  transfer  constant  (d1) 

•  Kd  -  Sorption  coefficient  (mg/L) 

•  Beta  -  First  order  dissolution  constant  (d1) 

•  Cs  -  NAPL  solubility  (mg/L) 


Figure  A.  9:  Chemical  Reaction  Package 


MT3D  =>  Run  RT3D 


Read  Results 


Breakthrough  Curve  Data 
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Display  =>  Plot  Wizard  -  Select  Time  Series  -  Next  -  Select  Show  All  -  Finish 

Once  chart  appears,  right  click  on  axis  and  select  Export/Print  -  Select  Text/Data  -  Select 
File  (Pick  destination  by  clicking  on  Browse)  -  Select  Export  -  Select  Maximum 
Precision  -  Select  Export 

Breakthrough  curve  data  can  be  generated  in  Excel  from  Exported  data  file. 

Mass  Data 


After  simulation  is  complete  Save.  It  is  important  data  is  gathered  from  the  desired  time. 


To  ensure  desired  time  is  selected,  select  species  in  home  screen,  then  time  under  Project 


Explorer  shown  in  Figure  A.  10. 
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Figure  A.  10:  Select  data  set  for  mass  analysis 
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MT3D  =>  Basic  Transport  Package  =>  Select  Species  >  Starting  Concentration  =>  3D 
Data  Set  ->  Grid  shown  in  Figure  A.  1 1 


Figure  A.  11:  Import  3D  data  set 

Once  clicked,  box  shown  below  will  appear.  Select  your  RT3D  data  set,  then  the  species 
you  want  to  study.  It  is  important  to  ensure  the  time  selected  is  the  end  time  in  the 
simulation  shown  in  Figure  A.  12. 
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Figure  A.  12:  Select  data  set  for  import 

Grid  will  then  be  filled  with  data.  Copy  pertinent  data  into  Excel.  Click  Cancel  after  data 
has  been  removed. 


Simulation  with  No  Source  Present 

The  following  instructions  are  for  after  the  source  has  been  removed.  Save  file  as  a 
different  name,  rerun  MODFLOW.  Delete  all  entries  in  Source/Sink  Mixing  Package. 
Follow  directions  above  for  Mass  Data,  however  do  not  click  cancel  once  data  sets  are 
imported  into  grid,  click  OK.  Change  Stress  Period  to  desired  length,  then  run  RT3D. 
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